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SUMMARY 


This is a final report on the tasks supported by NASA Langley Research Center un- 
der Grant NAG1-1021, Analysis, Preliminary Design and Simulation Systems for 
Control-Structure Interaction Problems. When the proposal was submitted, it was 
intedned to be a three-year program, with its first-year effort to be concentrated on soft- 
ware aspects of control-structure interaction(CSI) analysis. Due to the termination of the 
grant at the end of its first-year period, appropriate adjustments had to be made in order 
to make most out of the grant. The accomplishments from the one-year effort include: 1) 
tLe delivery of a research-level CSI analysis software that runs both on SUN Workstations 
as well as Alliant shared- memory parallel machines to Dr. W. K. Belvin of NASA/Langley 
Research Center; 2) three presentations at conferences, one in a bound book publication, 
the second in the 1990 AIAA Guidance and Control Conference Proceedings, and the 
third to appear in a NASA/DOD/JPL proceeding on system identifications. Two of these 
papers are being prepared for submittal to journal publications. 
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Summary 


This report summarizes research work on the implementation of a partitioned numerical 
algorithm for determining the dynamic response of coupled structure/controller/estimator 
finite-dimensional systems. The partitioned approach leads to a set of coupled first and 
second-order linear differential equations which are numerically integrated with extrapola- . 
tion and implicit step methods. The present software implementation, ACSIS, utilizes par- 
allel processing te chni ques at various levels to optimize performance on a shared-memory 
concurrent/vector processing system. The current work also generalizes the form of state 
estimation, whereby the Kalman filtering method is recast in a second-order differential 
equation equivalent to and possessing the same computational advantages of the structural 
equations. As part of the present implementation effort, a general procedure for the design 
of controller and filter gains is also implemented, which utilizes the vibration character- 
istics of the structure to be solved. Example problems are presented which demonstrate 
the versatility of the code and computational efficiency of the parallel methods is exam- 
ined through runtime results for these problems. A user’s guide to the ACSIS program, 
including descriptions of input formats for the structural finite element model data and 
control system definition, can be found in Appendix A. The procedures and algorithm 
scripts related to gain design using PRO-Matlab axe included in Appendix B. In Appendix 
C, a stability analysis for the partitioned algorithm is presented which extends previous 
analysis to include observer dynamics, leading to a clearly definable stability limit. The 
source code for the parallel implementation of ACSIS is listed in Appendix D. 


1.0 Introduction 

The present work on the implementation of a partitioned transient analysis algorithm for 
the simulation of linear Control-Structure Interaction (CSI) problems has concentrated 
on four major areas. The initial software implementation emphasized the user-friendly 
aspect and a structural dynamics-oriented interface for experienced practitioners of finite 
element analysis progr ams . Another reason for a new implementation of the algorithm 
was that the initial architecture of the CS3 software testbed developed by Belvin and Park 
[1-2] had little provision for effective parallelization. CS3 also included extensive links to 
optimization and optimal control algorithms which were not central to the current work 
and proved to be a further hinderance. This new software implementation was designated 
ACSIS, for Accelerated Control-Structure Interaction Simulation, and led to significant 
improvements in speed for particular problems on conventional serial processors due to 
simpler and more economical storage of primary variables. ACSIS is also versatile in its 
usage of a general Timoshenko beam element with pin release capability and shear correc- 
tion factor adjustment, and control system definition via a single data file. A preliminary 
Users Guide was developed for ACSIS, and the input formats were made compatible with 
pre/postprocessing software developed for Sun and Silicon Graphics computers so that 


future versions using parallel techniques via element mesh domain decompositions can rely 
on the same X- Window-based I/O utilities. Table 1 documents runtime comparisions of 
ACSIS and CS3 on a sample 48 DOF problem simulated on a Sun 3/260 workstation using 
a floating point accelerator. 


2.0 State Estimation via Second-Order Kalman Filtering 

One restriction of the CS3 software testbed was the form of dynamic observer equations 
used in the partitioned algorithm as developed in [2]. However, Belvin and Park [3] showed 
that a general Kalman filtering type of state estimation was not only possible in the parti- 
tioned solution, but could be implemented at a very small additional cost in computations 
by a slight modification in the way the dynamical equations of the plant are cast into 
first-order form for filter gain design. The methods employed in [3] have been successfully 
implemented into ACSIS, thereby enhancing the code’s ability to handle a wide variety 
of state estimation schemes. This is important as the application of optimal control tech- 
niques to the state estimator problem leads to this more general form, and as the restric- 
tions of the observer used in CS3 typically meant either discarding part of the observer 
gain parameters, resulting in a loss in system performance. 

ACSIS retains the option of using the more resricted observer form, as well as simple 
full state feedback (no dynamic compensation). It is clear from examining the respective 
equations, that for structures with stiffness- proportioned damping, the only additional 
calculations required for the K alman filter is the multiplication of the Li gain matrix by the 
predicted state estimation error vector j (see equations 25 of [3]). This is not a significant 
portion of the computations required at each integration step, as the dimension of 7 is small 
(number of sensors). By far the major costs are for computing an internal force of the 
form Kq, and backsolving the factored integration matrix for the estimated displacement 
states. This is verified numerically in the ACSIS results of Table 1 (using a restricted form 
of observer) and version (Al) of Table 2. The model used in Table 1 is roughly comparable 
to the 54 dof truss in Table 2, and both show the additional simulation time due to state 
estimation is roughly the same as an additional transient analysis. Finally, the Kalman 
filtering equations do not lead to any additional complications in parallel implementation of 
the overall algorithm as compared to the more restricted second-order observer developed 
for CS3. 
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3.0 Parallel Implementation of ACSIS 

A primary emphasis in our work dealt with the optimization of the software implemen- 
tation on a concurrent processing system. The platform chosen (primarily due to avail- 
ability, initially at CU and later at NASA LaRC) was the Alii ant FX/8 shared memory 
multiprocessor system with 8 parallel processing units and vectorization capabilities. Ver- 
sions of the software have been ported to the Alliant and compiled using the Concentrix 
FX/Fortran optimizing compiler, which has available options for automatic vectorization 
and concurrency of standard, problem-independent parallel computations. 

The partitioned CSI algorithm has three primary levels of parallelism in its numerics which 
can be exploited. At the highest level is the integration of the second-order dynamical plant 
(structure) and filter (estimator) equations, which as designed are of roughly equivalent 
size. Through the algorithm, these systems are effectively decoupled and independent 
at each discrete time step, and thus may be handled in parallel by invoking a compiler 
directive in the main program, which calls the respective subroutines simultaneously and 
handles re-synchronization of the execution upon their return. 

At a lower level of the algorithm, the structure and state estimator equations exhibit 
symmetric, second-order forms typical of linear structured dynamics problems. It is well 
known that computations related to the formation of the internal force vector, Kq, can be 
re-implemented at an element level [4], which, through decomposition of the element mesh 
[5], can be handled in parallel within exclusive subdomains. This technique becomes par- 
ticularly attractive for larger problems on more massively parallel systems; in the present 
work, the element-by-element (EBE) technique was not as effective as other types of paral- 
lelism. It should be noted here that the partitioned algorithm employs implicit integration 
methods, leading to systems of algebraic equations which are factored and solved using di- 
rect, rather than iterative, numerical methods. Therefore, formulation of the internal force 
vector is needed only in the formation of the known (right hand side) vector for integration 
of the plant and filter equations. The alternative to EBE computations is multiplication of 
the relevant displacement vector with the global stiffness matrix stored in profile (skyline) 
form. To “simulate” the advantages of local memory typical in large-scale parallel process- 
ing systems, the computed element stiffnesses were saved in shared memory for the EBE 
calculations, thus avoiding the need to recompute this data at each integration step. In 
addition, a low-overhead automatic element domain decomposition was provided for the 
parallel EBE method. 

The final and lowest level of parallelism is obviously that of the basic matrix computations 
such as addition, multiplication, etc. These numerical operations me inherent in nearly 
all areas of the software implementation, including problem preprocessing. With the very 
capable optimizing compiler available on the Alliant system, this parallelism was exploited 
through vectorization and concurrency of the nominal source code using the FX/Fortran 
compiler run-time options -0 -DAS -alt. The performance of the resultant executable 
code was examined using the A lli ant’s profiling capabilities, and changes to the nominal 
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source code, remaining compliant to F77, to maximize the identifiable concurrency and 
thus enhance the resultant speed. This was particularly useful in the profile matrix/ vector 
multiply operation, whose speed is critical to the overall program performance. 


4.0 Problem Descriptions 


Three structural dynamics problems were developed for code testing at various levels of 
complexity. All three problems have the following common features: simulations consisted 
of 1000 integration steps and employed a stiffness-proportional damping. The damping 
was not needed for algorithm stability, but to ensure consistency between the examples 
and because the existence of damping in the plant equations has a strong influence on 
program speed. The Kalman filter models were of second-order form [2] and equivalent in 
size to their respective plant models. For controlled simulations, the control system began 
operating after 100 integration steps, and all gain matrices were full (i.e. all model states 
influence all actuators and are influenced by all sensors). Additional specific information 
for the problems follow. 


A. Axial Vibration of Elastic Bar (Spring Model) 


# Nodes: 

# Elements: 
# Degrees of freedom: 

# Actuators: 
# Sensors: 

Disturbance: 


3 free, 1 fixed 
3 
3 
1 
1 

Initial displacement 


B. Planar Vibration of Space Truss (Truss Model) 


# Nodes: 

# Elements: 
# Degrees of freedom: 

# Actuators: 
# Sensors: 

Disturbance: 


18 free 

33 

54 

4 

6 

Bang-bang type sinusoidal applied force 


C. General 3D Vibration of EPS Satellite Reboost (EPS7 Model) 


# Nodes: 

# Elements: 
# Degrees of freedom: 

# Actuators: 
# Sensors: 

Disturbance: 


97 free 
256 
582 
18 
18 

Bang-bang type square wave applied force 
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As can be seen, the problem sizes are roughly three different orders of magnitude, with 
corresponding increases in the sizes of the control systems. Appendix B includes routines 
using Matlab for the design of control and filter gains which were used for the control 
system design of all three example problems. An illustration of the EPS model is shown 
in Figure 1. 


5.0 Performance Assessment 

Table 2 compares CPU runtimes on the Alii ant computer (using the UNIX "time” com- 
mand) for distinct versions of the software. Version (Al) is the nominal F77 program code 
compiled without any performance-enhancing options, while version (A2) invokes auto- 
matic vectorization and concurrency of low-level, problem-independent computations such 
as vector addition and inner products. The performance improvements are significant, 
especially for the large EPS7 model, where the speed-up factor is 35-37. 

Version (A3) also uses the compiler options from (A2), but in addition has a compiler 
directive added to the main program which allows the plant and filter integration subrou- 
tines to be called in parallel. This does not affect the transient response results as that 
analysis option bypasses the altered code, but for controlled response there is some effect 
on performance. If filtering is used, which results in a significant increase in computation, 
the directive can lead to some increased speed as can be seen for the spring and truss prob- 
lems. There can also, however, be a reduction in performance as compared to (A2) if the 
finite amount of processors and vector units axe used in a less efficient way. This appears to 
be the case for the large EPS problem, where the "overhead” introduced by the directive, 
and its effect on processor assignment, is greater than the improvement generated by the 
manually-invoked parallel construct. 

Table 3 shows CPU run times for ACSIS using E-B-E computations, which, as mentioned 
previously did not lead to better performance on the example problems using the Alliant 
system. This appears to be due to the lack of effective vectorization of the individual 
element computations when forming the internal force via the EBE method. To deter- 
mine whether the EBE method effectively lead slower speeds through increased numbers 
of computations, version (A2) (see above) was altered by removing compiler optimization 
of the profile matrix/vector multiply routine (the alternate method to EBE). The resultant 
runtimes matched almost exactly with those of version (A4), leading to the conclusion that 
both methods require roughly the same amount of computations, but differ in how they 
can be optimized on the Alliant system. The matrix/vector multiply operation, in this 
environment, can exploit both vectorization and concurrency through the compiler’s per- 
formance options; this can be examined in the compiler output. The EBE computations, 
at the element level, do not vectorize because the parts of the global displacement and force 
vectors being operated on per element are not contiguous. In version (A5), the elements 
within each subdomain of the mesh are computed in parallel, leading to some performance 
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Overall, versions (A2) and (A3) provide the best code performance for the hardware avail- 
able. Parallelizing the observer and structure (A3) leads to mixed results; improvement 
for the small spring and truss problems, but not for the large EPS model. Element-by- 
element computations do not improve code performance over compiler optimization via 
vectorization and concurrency for this platform. Reimplementation of the algorithm lead 
to a 5:1 improvement over the CS3 testbed software on a serial computer (Table 1). Fur- 
ther optimization of ACSIS on the Alliant FX/8 lead to an additional 30:1 improvement in 
runtimes for large-order systems such as the 582 dof EPS model. Time history responses 
of selected variables for the example problems axe shown in Figures 2 through 10. For 
Figures 8 through 10, the u x ,u y ,u z displacements plotted axe located at node 45, which 
is located at the vertex of the large antenna of the EPS model (see Figure 1). 


6.0 Conclusions 

The present work has demonstrated the efficiency of a streamlined simulation code for 
the analysis of large-order CSI systems and the viability of the continuous second-order 
Kalman filtering equations for state estimation. These methods show the versatility of 
the partitioned CSI integration algorithm and the promise of its application to real-time 
simulation. It is evident, however, that the use of element-by-element computational tech- 
niques requires the development of innovative algorithms for effective implementation on 
massively parallel processing systems. Future work in this area will include integrating 
algorithms for on-line system identification and applying these capabilities to the problem 
of real-time control. 
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Simulation 

CS3 

ACSIS 

Transient 

439.2 

98.8 

Full State 

Feedback 

688.2 

181.5 

FSFB with 

Observer 

1156.7 

282.3 


Table 1: Comparison of Rimtime Speeds for CS3 and ACSIS 

on a Sun 3/260 System 




(Al) 

(A2) 

(A3) 


Problem 

Nominal 

Compiler 

Parallel 

Model 

Type 

Code 

Optimized 

Observer 

3 DOF 

Transient 

6.6 

2.1 

2.1 

Spring 

FSFB 

8.0 

3.3 

3.3 


K. Filter 

12.3 

3.5 

3.3 

54 DOF 

Transient 

78.2 

5.7 

5.6 

Truss 

FSFB 

97.1 

9.4 

10.2 


K.Filter 

170.7 

13.0 

10.7 

582 DOF 

Transient 

3506. 

98.6 

100.3 

EPS7 

FSFB 

7040. 

190.2 

294.5 


K. Filter 

n/a 

284.2 

312.5 


\ 


Table 2: CPU Results for Versions of ACSIS 
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Model 

Problem 

Type 

(A4) 

E-B-E 

Computation 

(A5) 

Parallel 

E-B-E 

(A6) 

Parallel 
Obs. & EBE 

3 DOF 

Transient 

3.8 

3.3 

3.3 

Spring 

FSFB 

4.9 

4.4 

4.9 


K. Filter 

6.6 

5.6 

5.0 

54 DOF 

Transient 

31.7 

13.0 

13.0 

Truss 

FSFB 

35.5 

16.9 

35.6 


K.Filter 

62.6 

27.3 

36.2 

582 DOF 

Transient 

391.7 

153.9 

n/a 

EPS7 

FSFB 

485.9 

245.9 

n/a 


K. Filter 

n/a 

n/a 

n/a 


Table 3: CPU Results for ACSIS with EBE Computations 
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Figure 1: EPS Finite Element Model 
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Truss Model: Open Loop Transient Response 
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Figure 5: Truss Transient Response 
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Truss Models Full State Feedback Response 
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Figure 6: Truss FSFB Response 
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Truss Model: Controlled Response w/Kalman Filter 
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Figure 7: Truss Response w/Filter 
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Accelerated 

Control 

Structure 

Interaction 

Simulation 


Introduction 

ACSIS is an analysis program for full-order simulation of control-structure interaction 
(CSI) problems. The CSI simulation is carried out using a partitioned analysis pro- 
cedure which treats the structure (or plant), the observer, and the controller/observer 
interaction terms as separate entities. This procedure allows ACSIS to maintain rela- 
tively small, sparse matrix equations when compared to the process of assembling the 
computational elements into a single set of equations of motion and solving simultane- 
ously. Although this software can also carry out modal analysis, the transient response 
and CSI simulation are done in real space using the entire finite element model. (For 
more information, see Ref. 2 of report, p. 14) 

The ACSIS program is run using two previously prepared data input files and interactive 
input of run options. The two data input files are the finite element input file and the 
controller definition file. 
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Interactive Options 

The run options are as follows, note that no defaults exist, data must be input each 
time it is requested for each item requested except during a background run. Files may- 
be given any names, example names are given only to match the truss example at the 
back of this manual. 

Please input analysis type: 

This option is for selecting modal analysis, CSI simulation or the transient response of 
a structure. Not ail interactive inputs are required for each analysis type. 

Do you wish to save an input file? (y or n) 

This option creates a file which saves all the interactive input options. ACSIS can be 
background run with minor option changes by editing this file and then directing the 
screen input to this file. This file is very different for each analysis type . (To do this 
run with the example names after running acsis.exe once interactively, the command 
would be: acsis.exe <INP_truss> &. 

Name of save input file? (filename) 

This question asks for the name under which to save the interactive input, example 
name: INP -truss 

Finite Element Model Input File Name (filename) : 

This file should contain all the finite element nodes, mesh, materials, properties, lumped 
inertias, fixations, and initial velocity and displacement conditions. It must be prepared 
in advance in the card format specified later, example name: FEM_truss 

Number of modes desired? - 

This only appears in the modal analysis to request the number of modes to be output. 
The actual analysis is carried out with double this number of modes for accuracy. 

Controller Definition File Name: 

This only appears for the CSI simulation. The file should contain the actuator and 
sensor locations, control gains, and observer gains. It also must be prepared in advance 
in card format, example name: C0N_truss 

Please input type of control:: 

This option is for selecting the form of control law equations used in the CSI simulation. 
Full state feedback uses the current states of the plant (structure) to determine the 
control via constant gain matrices. Second-order observer uses only the L2 filter gain 
matrix of a Kalman filtering type of state estimation where the state variables axe 
position and velocity. The Kalman filter option allows the use of a full set of filter 
gains, but the gain design must come from a alternate variable casting using position 
and generalized momenta. 
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Initial time, final time, control-on time, step size: 

Data format is four col umns for CSI, but only three columns for transient response since 
the control-on time part is deleted. 

Forcing function ID, scale factor, damping coeff- a,b: 

In order for any time dependent forcing function to be easily implemented despite vari- 
able time step sizes, the forcing functions must be entered into the subroutine f orces.f . 
The format is to assign each new forcing function an ID number in a elseif statement. 
An example of forces.f is included at the back of the manual. Then forces.f must 
be recompiled into an object file and acsis.exa relinked. This is easily done by typing 
make which detects changes to the fortran files and does only the necessary compiling. 
A particular forcing function is chosen by entering its ID number in the first column, 
in addition the force can be scaled by a constant using the scaling factor in the second 
col umn - The Rayleigh damping coefficients a and b are the third and fourth columns. 

Phase lag fix? (y or n) : 

Specifies whether to include an extra iteration of control and sensor state prediction 
at each integration step to improve accuracy. Not generally needed unless the user is 
investigating the source of instability in a simulation and needs to test the sensitivity 
of the response to the partitioned algorithm’s extrapolation method. 

Gain scale factors (4 toted) : 

These scaling factors are, in order: the FI control gain matrix (displacement), F2 control 
gain matrix (velocity), LI and L2 state estimator filter gain matrix. This question 
appears only for CSI. 

ACSIS has two types of output options. The first is the displacement or velocity motion 
of up to twenty separate degrees of freedom. Interactive questions are, ‘Number of 
displacement results to output (max 10)’ followed by as many ‘Input node #, 
dof for displacement output #’ as necessary. Then these repeat for velocity results. 
The name of the file where the output will be stored is requested by ‘Output file 
name? (filename).’ example name: OUT.truss. The second output option saves 
the displacement of all nodes at any time step where output is sent (see next entry). 
The format is suitable for animation of the entire structure. Question asks ‘Animation 
Output? (y or n) ’ then for an animation filename if necessary. 

Send output every how many steps? 

This option affects both output options to reduce the size of the output and animation 
files. It causes output to only be sent after a integer number of time step iterations. To 
get output each time step, simply enter 1. 
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Finite Element Input File 

The finite element input file consists o£ title cards followed by columns of data- The title 
cards can be in any order, but they must be all capitalized with the appropriate number 
of columns of data for each card. The program reads rows of data until encountering 
a new card. Data which represents an integer value may be entered with a decimal 
point while real data may be entered without a decimal point as necessary. Any line 
beginning with a * anywhere in the file is ignored and can be used to insert comments. 
Any blank line will result in a read error. 

NODES 

Each node must be defined on a separate line. Columns can be separated by any number 
of spaces or a comma. Data format is four.columns: node number, x-coordinate, y- 
coordinate, z-coordinate. 

TOPOLOGY 

Each element must be defined on a separate line. Truss elements require two nodes then 
two columns of zeroes. (Truss elements would also require pin releases in ATTRIBUTES 
below.) Beam elements require two nodes then a third reference node representing a 
point in the xz plane of the beam and then a column of zeroes. Element type refers to 
finite elment formalation. (Currently only type 1 = timeshenko beam element is now, 
available.) Data format is six columns: element number, element type, node #1, node 
■#2, node #3, node #4. 

ATTRIBUTES 

Each element is characterized by an ID number from each of the MATERIAL and PROP- 
ERTIES cards below. Each element also has six pin release codes. The first code is for 
longitudinal stiffness, the second code is for torsional stiffness, the third and fifth codes 
are bending stiffness at each end in the y direction and must be the same value, and the 
fourth and sixth codes are for bending stiffness in the z direction and also must be the 
same. (0 is stiff, 1 is released.) Data format is nine columns: element number, material 
type, property type, and six pin release codes. 

MATERIAL 

The material data is formatted in four columns: material type, Young’s modulus, shear 
modulus and density of the material. 

PROPERTIES 

The properties data is formatted in six columns: property type, cross sectional area of 
the element, I y * I z , I y , I z , shear shape factor SSF2, shear shape factor SSF3. 

FIXITY 

Nodes with any fixations are defined here in the finite element file. The nodes must be 
entered with the fixity of all six of their DOF’s, restrained or not. (1 is restrained, 0 is 
free) Data format is seven columns: node number, x, y, z, 4> z , <p y , <j> z . 
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INERTIA 


All lumped inertias must be entered with each separate DOF on an individual line. 
Therefore a single node could take up to six lines to define. Data format is three 
columns: node number, DOF number (1-6), and value of inertia. 

INITIAL CONDITIONS 

All initial displacement and velocity conditions are entered into the finite element file. 
Data format is four col umns : node number, DOF number (1-6), initial displacement, 
and initial velocity. 

END 

End of file. 

Controller Definition Cards 

The controller defintion file also consists of title cards followed by data entry. Each card 
must be followed by the appropriate number of columns of data and in some cases the 
appropriate number of rows. Integers can be entered in real format and vice versa if 
necessary. Any blank line will result in a read error. 

NACT 

Number of actuators in the entire control system. 

BMAT 

This entry creates the actuator position matrix or B matrix. There should be one 
row for each actuator, data format is four columns: node number, DOF number (1-6), 
actuator number, and sensitivity. 

NSEN 

i 

Number of sensors in the entire controls system. 

HDMA 

This entry creates the matrix of displacement sensor locations. One row per displace- 
ment sensor, data format is four columns: node number, DOF number (1-6), sensor 
number, and sensitivity. 

HVMA 

This entry creates the matrix of velocity sensor locations. One row per velocity sensor, 
data format is four col umns : node number, DOF number (1-6), sensor number, and 
sensitivity. 

F1GA 

This is a list of the Fl or displacement control gains. The data format is four columns: 
node number, DOF number (1-6), actuator number, and value of gain. 
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F2GA 


This is a list of the F2 or velocity control gains. The data format is four columns: node 
number, DOF number (1-6), actuator number, and value of gain. 

L1GA , 

This is a list of the state estimator Ll filter gains. The data format is four columns: 
node number, DOF number (1-6), actuator number, and value of gain. 

L2GA 

This is a list of the state estimator L2 filter gains. The data format is four columns: 
node number, DOF number (1-6), actuator number, and value of gain. 

END 

End of file. 

Examples 

This section includes all the files and procedures needed to run all three analysis on a 
simple elastic bar problem. The naming of the files is a simple and easy to remember 
system, however no particular format is necessary. 

The finite element file was created simply by typing the node locations, connectivity 
(topology), etc. with a text editor. 

File: FEM_spring 


MESH 

* 

NODES 

1 0.00 0.00 0.00 

2 1.00 0.00 0.00 

3 2.00 0.00 0.00 

4 3.00 0.00 0.00 

TOPOLOGY 

1 112 0 0 

2 1 2 3 0 0 

3 1 3 4 0 0 

ATTRIBUTES 

1 11011111 
2 11011111 
3 1 1011111 

MATERIAL 

1 1000 . 0.0 0.0 

PROPERTIES 
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1 1.00 0.00 0.00 0.00 0.0 0.0 

FIXITY 

1 111111 
2 0 11111 

3 0 11111 

4 0 11111 
INERTIA 

2 1 0.100 

3 1 0.100 

4 1 0.100 

INITIAL 

3 1 0.100 0.000 

END 


The modal analysis only requires the number of modes desired in addition to the finite 
element file. The file INP-springO documents the interactive inputs used in the modal 
analysis. The results axe saved in file EIG .spring. 


File: INP_spring0 


-1 

n 

* ACSIS input file, two lines above are 

* analysis type and save input file. Do 

* not change then by editing this file. 

* Finite element input file?(f ilenane) 
FEM_spring 

* Number of inodes desired? 

3 

* Output file? (filename) 

EIG.spring 


File: EIG_spring 


SUBSPACE ITERATION ROUTINE 

NB OF EIGENVALUES =» 3 

NB OF VECTOR* 3 

NB OF D0F» 3 
TOLERANCE- 1.000E-04 

NB OF RIGID MODES- 0 


ITERATION NO 1 
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l.OOOE+OO 1 . OOOE+OO l.OOOE+OO 

ITERATION NO 2 
1 .297E-14 3 . 194E-14 3.899E-14 

EIGEN ANALYSIS RESULTS: t < 

RADIAL - CYCLIC 

MODE EIGENVALUE FREQUENCY FREQUENCY 


1 

1980.6 

44.504 

7.0831 



2 

15550. 

124.70 

19.846 



3 

32470. 

180.19 

28.679 



EIGENVECTORS : 





1 

2.3305 

1.8689 

-1.0372 

0. 

0 . 

2 

1.8689 

-1.0372 

2.3305 

0 . 

0 . 

3 

1.0372 

-2.3305 

-1.8689 

0 . 

0 . 


MASS MATRIX DIAGONAL: 

213 l.OOOOOOOOOOOOOD-Ol 
312 l.OOOOOOOOOOOOOD-Ol 
411 l.OOOOOOOOOOOOOD-Ol 


To run the transient response of the structure, a forcing function or initial condition 
would be needed to excite the structure. An initial condition would be added to the 
finite element file. A forcing function must be added to f orces.f with a new ID number, 
then this number given as interactive input. In this case, an initial displacement acts 
on the second degree of freedom, which is defined in the finite element model input file. 
The file INP jspringl documents the interactive inputs used in the transient analysis. 

File: INP_springl 


-3 

n 

* ACSIS input file, two lines above are 

* analysis type and save input file. Do 

* not change them by editing this file. 

* Finite element input file? (filename) 
FEM. spring 

* Initial, final, step size? 
0.00000000 1.00000000 0.00100000 

* Forcing function, scale f, damping a,b? 
0 0.000000 0.00000000 0.00002000 

* Output file name? (filename) 

0UT_spring 

* Number of displacement outputs? 

1 

3 1 

* Number of velocity outputs? 
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0 

* Send output every how many steps? 

1 

* Send animation output?(y or n) 
n 

In addition to the finite element file, the interactive input, and an excitation, CSI 
simulation requires a controller definition file. A full state feedback controller for the 
truss structure is defined in COS_spring. 


File: COH_spring 


NACT 

1 

BMAT 

3 1 1 1.0 
F1GAIN 

211 193.31415000000 

311 1574.6315000000 

411 -1669.0460000000 

F2GAIN 

211 20.774256000000 

311 27.790403000000 

411 10.780212000000 

NSEN 
1 

HVMAT 
3 1 1 1.0 
L1GAIN 

211 8 . 7928909000000D-02 

311 -4. 880790 1000000D-02 

411 3 . 0102735000000D-03 

L2GAIN 

211 1 . 0380932000000D-0 1 

311 6.1410130000000 

411 -3.0608725000000 

END 



Then if acsis.exe is run interactively and the input is saved, the file INP_spring2 can 
be produced. The file could be edited to change the input files, the output file, the 
forcing function ED, the length of simulation, control-on time, etc. Then to run it again, 
type acsis.exe<INP_spring2> scr A. The scr is a scratch file which will store the 
screen output. 


File: INP_spring2 


2 
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* ACSIS input file, two lines above are 

* analysis type and save input file. Do 

* not change then by editing this file. 

* Finite element input file? (filename) 

FEM_ spring 

* Controller file name? (filename) 

CON. spring 

* Please input type of control: 

3 

* Initial, final, control-on, step size? 

0.00000000 1.00000000 0.10000000 0.00100000 

* Forcing function, scale f, damping a,b? 

0 0.000000 0.00000000 0.00002000 

* Phase lag fix?(y or n) 

n , 

* Gain scale factors (4 total)? 

1.00000000 1.00000000 1.00000000 1.00000000 

* Output file name? (filename) 

0UT_ spring 

* Number of displacement outputs? 

1 

3 1 

* Number of velocity outputs? 

0 

* Send output every how many steps? 

1 

* Send animation output? (y or n) 
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Introduction 


In order to retain simplicity in the ACSXS code for parallel implementation purposes, 
no control system design algorithms were included. Instead, using modal data out- 
put from the eigenmode analysis module of ACSIS, a procedure was developed using 
the Pra-Matlab and its Control System Toolbox, which includes algorithm scripts for 
optimal control solutions via the solution of an algebraic Ricatti equation. In order to 
accomodate large order dyn ami cal systems, the design is accomplished in the uncoupled 
normal modes domain, using the available lowest eigenmodes from ACSIS. 

The procedure begins by copying and editing the mode data output from ACSIS (see 
the listing for EIG_spring in Appendix A) into readable variables for input to Matlab. 
The typical approach was to create one file as a Matlab script (i.e. a ” .m” file), with 
the eigenvalues, eigenvectors, and mass/dof data from ACSIS at the beginning, followed 
by actuator and sensor influence matrices (related to the physical degrees of freedom), 
the objective function weighting matrices, and the function which calls other Matlab 
scripts to determine the solution. An example of the above input for the spring problem 
described in Appendix A is in file EIG_spring.m, which is listed below. Compare this 
the the ACSIS mode data output shown in Appendix A to see how the editing was 
accomplished, and the additional control design data added 


File: EIG_spring . m 


lam ■ [ 

1980.6 

15550. 

32470. 

]; 




vi®[l 

2.3305 

1.8689 -1.0372 

0. 

0. 


2 

1.8689 

-1.0372' 2.3305 

0. 

0. 


3 

m*[ 2 

1.0372 
1 3 

-2.3305 -1.8689 

1 .00000000000003-01 

0. 

0. 

3; 


312 1.0000000000000D-01 

411 1.0000000000000D-01] ; 

t»[vl(:,2:4)3; 
qb«*zeros(3,l) ; 
qb (2,1) *1.0; 
hd=zeros(l ,3) ; 
h.v»zeros(l,3) ; 
hv(l,2)®1.0; 
qv»[l; .5; .1] ; 
q=diag( [qv ; qv] ) ; 
r*.0001*eye(l) ; 

[f l,f2]=alqr(lam,t,a,qb,q,r,3) ; 
qv*Cl;l;l] ; 
q*diag([qv;qv3) ; 
r*100*eye(l) ; 

Ckl,k23*akf (laa,t,a,hd,bv,q,r,3) ; 
save flout fl /ascii 
save f2out f2 /ascii 
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save klout kl /ascii 
save k2out k2 /ascii 


The scripts mlqr.m and mkf .m were written to accept as input the vector of eigenvalues, 
the eigenvector matrix (orthogonal vectors stored in columns), a matrix of node, com- 
ponent, d.o.f, mass data, and the actuator (or sensor) influence matrix and weighting 
matrices. The scripts output the gain results in four- column arrays, with one gain per 
row, and the corresonding node number, displacement component, and actuator (or 
sensor) identification. The top-level problem script (listed above) then saves the output 
in external files and the analysis is complete. The design scripts (listed below) also 
include checks on controllability and observability of the system based on the modal 
data and influence matrices defined, and produce plots of the closed loop poles resulting 
from the gain design. This aids the analyst in assessing the expected performance (and 
stability) of the exact system before moving the data back to ACSIS for simulation. 
The script contrank.m finds the rank of the controllability matrix through iterative 
rank calculations of submatrices so as to avoid the iilconditioning experienced in the 
full matrix. This is botfr.faster and more accurate for determining whether a particular 
actuator placement has full control of the included structural modes. 

File: mlqr.m 


f unct ion [? lout , F2out] =mlqr (lam , t , m , qb fl , R , nmode ) 

7 . 

7 . Controller gain design for second-order 
X structural system via given eigenmodes. 

X Gains are transformed to be coefficents 
X of structural variables (disp .velocity) ; 

X i.e. plant is second-order, of size ndof. 

7 . 

V, Arguments : 

7 . 

X lam: vector of eigenvalues (nmode z 1) 

X t: matrix of eigenvctors (ndof x nmode) 

X m: mass diagonal and dof mapping info (ndof x 4) 

X qb: actuator position influence matrix (nact x ndof) 

X Q: optimal design state weighting matrix (2*nmode x 2*nmode) 
X R: optimal design feedbk weight, matrix (nact x nact) 

X 

X Flout: FI gain matrix for 2nd-order plant 

X F2out: F2 gain matrix for 2nd-order plant 

7 . 

X Written by K.F. Alvin 
X 

format short e 
nmodmax“length(lam) ; 

[ndof , nact] s s ize (qb ) ; 

X 

X Variables : 

X mass: mass matrix 
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% A: state transition matrix 

% B: actuator influence matrix 

% G: control gain matrix 

X 

massd(m(: ,3))»m(: ,4) ; 
mass a diag(massd) ; 

A=» [zeros (nmode) .eye(nmode) ;-l*diag(lam(l :nmode)) .zeros (nmode)] ; 

Amax® [zeros (nmodmax) .eye(amodmax) ;-l*diag(lam) .zeros (nmodmax)] ; 

B* [zeros (nmode ,nact) ; t ( : , 1 : anode) 1 *qb] ; 

Bmax» [zeros (naodmax ,nact ) ; t ’ *qb] ; 

disp( ’Number of structural modes and actuators used: 1 ) 
disp ( [nmode ,nact] ) 

dispC’Rank of the controllability matrix: *) 
dispCcontrank(A.B)) 

disp (’Determining controller gains for given system. . . *) 

G«lqr(A,B.q.R); 

vmax a l.l*max(sqrt(laa)); 

axis ( [-smax , smax , -umax , smax] ) ; 

plot(eig(A-B a G) , ’♦*) 

grid 

title( ’Roots of controlled system*) 

hold 

pause 

x 

X partition gain matrix 

X 

G1 S G( : , i : nmode) ; 

G2=G( : , nmode* 1 :nmode+nmode) ; 

y. 

X Transform resultant gain matrices for use in 
X partitioned csi algorithm using second-order 
X structure (plant) equations. 

X 

disp (’Mapping ga:'ns bach to physical domain...’) 
fl a Gl*t(: ,l:nmode) ’♦mass; 
f 2=G2*t ( : , 1 : nmode) ’ ♦mass ; 

X 

X find modal damping ratios of controller for calulated gains: 

X 

Gmax»[f l^t ,f2*t] ; 
lambda a eig(Amax~Bmax*Gmax) ; 
plot (lambda , ’ ♦ ’ ) 
pause 

nfreq a sqrt(imag(lambda) ."2 + reed, (lambda) . *2) ; 
mdamp<*-real (lambda) ./nfreq; 

disp (’Resultant modal damping ratios for controller:’) 
disp([’ Damping *,* Damped Freq (rad/s) ’]) 
disp( [mdamp. nfreq. ♦sqrt (1-mdamp. ~2) .lambda]) 
bdamp a max(-2^mdamp . /nfreq) ; 

disp( ’Estimated mi ninn™ stiffness damping coefficient necessary’) 
dispC to stabilize residual modes due to gain roundoff accumulation: *) 
disp(bdamp) 

disp( ’Writing gains in node correspondence output form...’) 
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Flout a zeros(nact*ndof ,4) ; 

F2out=zeros (nact*ndof , 4 ) ; 
for i«l:nact; 
kmin 3 (i-l)*ndof+l; 
taaax 3 i*ndof ; 

Flout (taain:tanax, :) a [m(: ,1:2) ,i*ones(ndof ,1) ,f l(i,m(: ,3)) ’] ; 
F2out(taain:taaax, :) a [m(: ,1:2) ,i*ones(ndof ,1) ,f2(i,m(: ,3)) '] ; 
end; 

disp ( * Finished mlqr * ) 


File: mkf .m 


function IX lout ,L2out]*mkf (lam, t,m,hd,hv,Q,R, nmode) 

X 

X Kalman filter design for second-order 
X structural system via given eigenmodes 
X and transformed to independent 
X displacement/gen . momentum variable 
X casting for partitioned csi transient 
X analysis. See Belvin/Park paper for 
X filter variable definitions. 

X 

X Arguments : 

X 

X lam: vector of eigenvalues (nmode s 1) 

X t: matris of eigenvctors (ndof x nmode) 

X m: mass diagonal and. dof mapping info (ndof x 4) 

X hd: sensor position influence matrix (nsen x ndof) 

X hv: velocity position influence matrix (nsen x ndof) 

X Q: optimal design state weighting matrix (2*nmode x 2*nmode) 
X R: optimal design feedbk weight, matrix (nsen x nsen) 

X 

X Llout: LI gain matrix for 2nd-order filter 

X L2out: L2 gain matrix for 2nd-order filter 

X 

X Written by K.F. Alvin 

X 

format short e 
nmodmax a length(lam) ; 

[nsen, ndof] =size(hd) ; 

X 

X Variables : 

X mass: mass matrix 

X A: state transition matrix 

X G: noise influence matrix 

X C: output influence matrix 

X K: filter gain matrix > 

X 

massd(m(: ,3))=m(: ,4) ; 
mass=diag(massd) ; 

A 3 [zeros (nmode) ,eye(nmode) ;-l*diag(lam(l:nmode)) , zeros (nmode)] ; 
Amax 3 [zeros (nmodmax) , eye(nmodmax) ;-l*diag(lam) .zeros (nmodmax)] ; 
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G=eye(nmode+nmode) ; 

C*Chd*t(: ,l:nmode) ,hv*t(: »l:nmode)3 ; 

Cmax a [hd*t ,hv*t] ; 

disp (’Humber of structural modes and sensors used:*) 
disp ( Cnmode ,nsen] ) 

disp ('Rank of the observability matrix: ’) 
disp (eontr ank (A * , C * ) ) 

disp (’Determining filter gains for given system...*) 
K a lqeU,G,C,Q.R); 

X 

X partition gain matrix 
X 

Kl a K(l:nmode, :) ; 

K2=K (nmode+ 1 : nmode+nmode , ; 

Kmax a [K 1 ; zeros (nmodmax-nmode ,nsen) ; K2;zero3 (nmodmax-nmode , ns en)] ; 

X 

X find modal damping ratios of filter for calulated gains: 

X 

lambda a eig(Amax-Kmax*Cmax) ; 

plot (lambda, ’+*) 

hold 

pause 

al 3 ~real (lambda) ; 
bl a imag (lambda) ; 

disp (’Result ant modal damping ratios for filter:’) 
disp(C' Damping ’,’ Freq (rad/s) ’]) 
disp([al./sqrt(al.~2+bl.“2) ,bl]) 

X 

X Transform resultant gain matrices for use in 
X partitioned csi algorithm using second-order 
X Kalman filter approach. 

X 

disp( ’Mapping gains back to. physical domain...’) 
ll a t(: ,l:nmode)*Ki; 

12=mass*t ( : , 1 :nmode) *K2 ; 

disp( ’Writing gains in node correspondence output form;..’) 
Llout a zeros(nsen*ndof ,4) ; 

L2out a zeros(nsen*ndof ,4) ; 
for i a l:nsen; 
kmin a (i-l)*ndof+l; 
kmax a i*ndof ; 

Llout (fcmin : kmax , :) a Cm(: ,1:2) , i*ones (ndof ,1) ,ll(m(: ,3) ,i)] ; 

L2out (kmin : kmax , :)*Cm(: ,1:2) ,i*ones(ndof ,1) ,12(m(: ,3) ,i)] ; 
end; 

disp (’Finished mkf ’) 


File: eontr ank. m 


function maxrankacontrankCa.b) 
maxrank=0 ; 

[instate, nact} a size(b) ; 
i=0; 
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»* 


mat°b ; 

nevr ank=*rani (mat) ; 
while nesrank > maxrank 
maxrank»newrank ; 
i a i+l; 

mat* [mat , (a"i) *b] ; 
nevrank«rank(mat) ; 
if newrank a *nstate 
maxrank=newr ank ; 
and 

end 


Unfortunately, the external files created from Matlab with the gain results are written 
completely in terms of real numbers, whereas the first three columns are actually to be 
read by ACSIS as integers (they are used as indices). A separate utility was written to 
convert the format of these files; the source code is listed below. On Unix systems, the 
user simply assigns the standard input to be the current data file created by Matlab, 
and gives another file name for the standard output. The code is basically just a filter to 
chahge the three columns of indices to integers. The output cm then be pasted directly 
into the control definition file used by ACSIS. 


File: convcont .f 


program convcont 

parameter (NMAX a 100000) 
real*8 f(NMAX),v(4) 

integer node(NMAX) ,dof (NMAX) , act(NMAX) 
n=0 

100 read (*,*,err a, 200,end a 200) (v(i) ,i a l,4) 
n«n+l 

node (n) =>int (v ( 1 ) ) 
dof(n) a int(v(2)) 
act(n) =int(v(3)) 
f(n) »v(4) 
goto 100 

200 do 300 i»l,n 

print *, node(i) ,dof (i) ,act(i) ,f (i)' 
300 continue 

end 
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APPENDIX C 


Stability Analysis of a 
CSI Partitioned Simulation 
Algorithm with State Estimator 



The equation of the open-loop plant without passive damping in modal second-order 
form is , 

q + u 2 q = u (1) 

The controller uses a second-order observer to estimate the plant state, along with a 
full-state feedback control gain design. 

u = - ( rju 2 p + (up) 

p + u 2 p = u + ( 7 ^ 

7 -z-p 
z = q 


where q and p are the plant and estimator states, respectively, u is the control force, 7 
is the state estimation error, 2 is the sensor output, and 77, (,( are gain coefficients for 
position and velocity feedback, and the estimator filter. 

The partitioned analysis procedure uses a stabilized form of the control law and esti- 
mation error determination to reduce inaccuracies associated with the extrapolation of 
variables in the controller force prediction. A first-order filtering is achieved by taking 
the time derivative of (2a, c), 

u = — rju 2 p — (up 
j = 2 -p 


( 3 ) 


and then embedding the equations of motion through substitution for p. This leads 
to the following two'coupled, .first-order differential equations for the prediction of the 
control force u and state error 7. 





TJU 4 

0 


} p 


(*) 


Time discretization of ( 4 ) using an implicit midpoint rule leads to the following coupled 
difference equation: 


1 +• S(u 

6 


6((u 

1 + 6 ( 




J S(u 3 — T]U 2 
{ Su 2 


K‘-{ c r}p- 


( 5 ) 


where 6 = half-step size = j. Solving this equation requires knowledge of the plant (to 
obtain sensor output) and observer states. These values are extrapolated as: 


Using these equations to obtain u n+ » and 7 n+ j allows the plant and observer equations 
to be solved independently. Midpoint time integration of (1) and (2b) leads to the 
following equations: 

(1 + *V) <Z n+ * = + q n + Sq n 

(1 + SW) J> n+ * = 6 2 u n+ * +p n + Sp n + 6 2 £ 7 n+ * 

i n+ * = J («” +} - 5") 

P" +i -f(p" +i -P") ( 7 ) 

q n+1 = 2 q n+ * ~ q n 

p n+l = 2 q n+ $ - q n 

4 n+1 = 2 q n+ $ - q n 

p n+x = 2q n+ ± - q n 


Computational stability of the modal form of the CSI partitioned equations of motion 
using the aforementioned time discretization can be assessed by seeking a nontrivial 
solution of 



such that 

|A|<1 


(9) 


for stability. Subtituting (8) into (5-7), we obtain 


Jx = 0 

where 


Xi 

II 

ST 

+ 

u n+ i 

p n+ j 

p n+ % 

p" ?"} 

X 2 

= [*? +i 

7 n +i 

q n +i 

q n +» 

?” ?”] 


Ju J12 
_ J 21 ^22 


( 10 ) 


( 11 ) 


40 


10 0 

(rju 2 — S(u> 3 ) (1 + 6 £uj) 0 


0 -1 

0 0 


Jn = 


J 22 = 


0 



- 6 2 (1 + 6 2 u 2 ) 

0- 

- 1 

-6 

0 



0 

1 

“7 




1 

i 

7 

0 



0 

-2 


A + 1 

0 

0 



0 

0 



- 

-2 

0 




‘0 

0 0 0 

0 

O' 







0 

6(& 0 0 

0 

0 





J 12 = 


0 

—6 2 Z 0 0 

0 

0 






0 

O 

O 

KKJ> 

1 

0 

0 







0 

0 0 0 

0 

0 







.0 

0 0 0 

0 

0. 






- 

0 

0 0 0 

0 

0* 






o 

o 

<0 

1 

0 

1 



■w 



0 

o 

o 

n 

1 

0 

0 



J 

21 = 


0 

0 0 0 

0 

0 






0 

0 0 0 

0 

0 





• 

0 

0 0 0 

0 

0. 



' 1 

0 


0 

0 


0 


-1 ' 

-1 

(1 + so 

0 

0 


0 


0 

0 

0 


(1 + S 2 U 2 ) 

0 


- 

1 

-<S 

0 

0 


1 

7 

1 


1 

6 


0 

0 

0 


-2 

0 

A + l 

0 

. 0 

' 0 


0 

_2 


0 


A + 1 . 


-s 

C u 

0 

A + 1 


( 12 ) 


(13) 


(14) 


(15) 


A nontrivial solution to (10) is found from 

det J = 0 (16) 

which leads to the characteristic equation 

((1 - <5£) (1 + <5 3 Cw 3 - 6 2 tju 2 ) + <5f (1 + 8 2 u 2 ))z i 

+ (1 — 6£) -J- 6£ (l + S 3 ^u 3 — S 2 t 7u/ 2 ))z 3 

((1 — 6£) ( S 2 u> 2 + 6 2 tju> 2 ) + 6 2 u 2 (l + 6 3 £uj 3 — 8 2 rju 2 ) 

+6£ (6£u> + 8 2 u> 2 + 8 4 u 4 ))z 2 
(S 2 u> 2 (8 (u>) -f* 5£ ( 6 2 uj 2 4- S 2 rju> 2 ))z 

+S 2 u> 2 ( S 2 uj 2 + S 2 rju 2 ) = 0 


where 
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A = 


1 + z 
1 -** 


|A |<1 


Re(z ) < 0 


( 18 ) 


Thus, a test of the polynomial equation for possible positive real roots by the Routh- 
Hurwitz criterion indicates that the partitioned. approach as applied to the modal equa- 
tions give a computationally stable solution for no velocity feedback C = 0 provided 


h < 


y/fjU 


( 19 ) 




(20) 
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APPENDIX D 


ACSIS Source Code 

tf 
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/ 


File: Makefile 


.SUFFIXES : .f 

FFLAGS « 

.f .0: 

fortran 

.0 

-c $( FFLAGS) $*.± 



OBJS 3 

aesis.o 

aesisont.o 

addstf .0 

\ 


bcaa 3 d . 0 

forces .0 

input . 0 

\ 


pa&vmul . 0 

prepf ea.o 

profile.o 

\ 


read.o 

solver. 0 

aephlag.o 

\ 


zerovect.o 

Xu.o 

prepcon . 0 

\ 


control. 0 

secorder.o 

measure. 0 

\ 


eigens . 0 

singeig.o 

animout . 0 

\ 


stif f rc . 0 
kf liter. 0 

eatifvm.o 

renum.o 

\ 

acsis.«x«: 

S(OBJS) 



1 

fortran 

-0 acsis.exe 

$(FFLAGS) $(OBJS) 



File: shared. inc 


t 



c 

c — 

c 

C shared. inc (ACSIS database) 

C 

C — 

c 

c 

C Argument definitions: 


C 

ad amp: 

Rayleigh damping coefficient alpha 

C 

b: 

actuator location matrix (packed storage) 

c 

brow : . 

row number of corresponding real value in b 

c 

bcol: 

column number of corresonding real value in b 

c 

bval: 

number of nonzero values in b 

c 

bdamp: 

Rayleigh damping coefficient beta 

c 

confile: 

controller input file name 

c 

contype : 

id for type of control 

c 

coxy z : 

array of the x,y, and z components of each node 

c 

delta: 

one-half time step 

c 

delaq: 

one-half time step squared 

c 

ec: 

control prediction integration coefficient matrix 

c 

•mat: 

array of element material types 

c 

eo: 

observer construct matrix S in vector form 

c 

eprop : 

array of element property types 

c 

es: 

structure construct matrix S (H+delta«D+delsq*K) 

c 

etype: 

array of element types 

c 

elnum: 

array of element numbers for domain decomposition 

c 

f : 

vector of applied forces 

c 

fl: 

control gain matrix 
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c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 


f2 : 

t emf ile: 
forceid: 
gamma: 

S c: 

gk: 

go-' 

S» 5 

h: 

id: 

hdrow: 

hdcol: 

hdval: 

hv: 

hvroo: 

hvcol: 

hml: 

id: 

ix: 

inertia.: 

jdiag: 

11 : 

12 : 

mask: 

maaa: 

mat: 

mien: 

nact: 

ncsi: 

ndisont : 

ndof : 

ndomain : 

nal: 

neld: 

nap: 

naan: 

nvelout : 

nolag: 

outfile: 

out label: 

pin: 

pivot : 

prop: 

<15 

q.0: 

qalpha: 
qalphao : 
qbata: 
qbetao : 
qdot: 
qdotO : 
qe: 

qadot : 
r: 

scalef : 

stiff: 

tO: 

tc 

ti: 


control gain matrix 

finite element input file name 

identification number of forcing function 

atate correction force 

RHS vector for control prediction module 

SSS vector for Kalman filter momentum eqn 

RHS vector for obaerver module 

RHS vector for atructure module 

time atop size 

displacement sensor location matrix (packed storage) 

ros number of corresponding real value in hd 

column number of corresonding reed, value in hd 

number of nonzero values in hd 

velocity sensor location matrix (packed storage) 

rov number of corresponding real value in hv 

column number of corresonding real value in bv 

number of nonzero values in bv 

OOF slapping array: id(comp,node #)=Global DOF 9 

array of element connectivity and orientation 

array of concentrated inertias or lumped masses 

array of diagonal element addresses 

State estimator filter gain matrix 

State estimator filter gain matrix 

mass matrix N in reduced vector form 
array of different nmterials 

length of global matrices in profile vector storage 

actual number of actuators 

actual number of actuators and sensors 

number of displacement results to output 

actual number of degrees of freedom 

actual number of element domains (for dom. decomp.) 

actual number of elements 

array of number of elements in- each domain 

actual number of nodes 

actual number of sensors 

number of velocity results to output 

logical flag to signal corrector loop in measurement 

output file name 

array of output data requested 

array of element pin release codes 

Column pivoting info from FACTA 

array of different properties 

generalized displacement vector 

initial displacement condition 

gain scale factor for fi 

gain scale factor for 11 

gain scale factor for f2 

gain scale factor for 12 

velocity vector 

initial velocity condition 

state estimator displacement vector 

state estimator velocity vector 

Solution vector of control module {u > gamma}’ 

Scaling factor for forcing function 
stiffness matrix K in reduced vector form 
initial time 
control-on time 
final time 
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vector of control forces 


C 

C 

c 

c 

c 

c 

c 

c 

e 

e 

c 

c 


a: 


Parameter definitions 


MIX ACT: 
MAX CSX : 
MAID AT: 
HAXDOF: 
HAXBQM: 
HAXELE: 
HAXHLEH: 
HAXHODE: 


max. * of actuators 

max. combined 8 of actuators and sensors 
max. * of a&terials and properties 
aax. t of degrees of freed©® 
max. S of decomposition domains 
max. * of elements 

max. length of global vectors in reduced form 
sax. # of nodes 


parameter(HAXDOF=3000 , HAXACTsSO, HAXCSI=100) 
parameter(HAXH0DE=1000, HAXELE=3000, MAID 11=100) 
parameter (HAXHL£N=200000 , HAXD0H»50, HAXHZV=200) 

real* 8 to , tf , t c , h .delta , delsq , qalpha , qbeta , qalphao , qbetao 
real*8 q(HAXDOF) .qdot(HAXSOF) .qe(HAXDOF) .qedot(HAXDOF) 
real*8 n( MAX ACT) ,gamma(HAXACT) ,f (HAXDOF) ,r(KAXCSI) 
real*8 es(MAXMLEH) , eo(HAXHLEN) . ec(HAXCSI .HAXCSI) 
real*8 gs (HAXDOF) .go(MAXDOF) .gc(HAXCSI) .scalef ,pe (HAXDOF) 
real*? mass(HAXHLEH) ,stif (HAXHLEH) , ad amp, bdamp,gk( HAXDOF) 
r eal*8 eoxyz ( 3 , HAXHODE) .mat ( 6 , HAXD AT ) , prop ( 1 0 , MAXD AT) 
read *8 qO ( 6 , HAXHODE) , qdo 1 0 ( 6 , HAXHODE ) , inert ia ( 6 , HAXNQDE ) 
real«8 b(HAXHZV) .hd(HAXHZV) .hv(HAXHZV) ,estifm( 78.500) 
real*8 fi(HAXACT. HAXDOF) ,f2(HAXACT, HAXDOF) 
real*8 ll(HAXDOF.HAXACT) , 12 (HAXDOF .MAXACT) 
integer etype(HAXELE) , ix( 4, HAXELE) ,emat(MAXELE) .forceid 
integer eprop(HAXELE) ,pin(6 .HAXELE) ,id(6 .HAXHODE) 
integer mask(HAXHODE) .contype .brov(HAXHZV) .bcol(MAXHZV) 
integer hdrow(HAXHZV) .hdcol(HAXHZV) .hvrow(MAXHZV) 
integer hvcol(HAXHZV) ,bval,hdval,hvval 

integer ndof ,nact,nsen, ncsi, mien, jdiag(HAXDOF) .nnp.nel.neig 
integer neld(HAXDOH) .elnum(HAXELE.HAXDOH) , eldom(MAXELE) .ndomain 
integer outlabel(40) .ndisout .nvelout .pivot (MAXCSI) 
integer iadjcy (HAXHLEH) ,icount(HAXHODE+l) .perm(MAXHODE) 
integer xls (HAXHODE) 
logical animate, nolag 

character*32 f emfile.conf ile.outf ile.animf ile 


COHHOH 

CQHHOH 

COHHOH 

COHHOH 

COHHOH 

COHHOH 

COHHOH 

CQHHOH 

COHHOH 


COHHOH 


/FILES/ f emf ile , conf ile ,outf ile , outlabel , ndisout , 
nvelout , animf ile , animat e , no lag 
/TIHERS/ tO, tf.tc.h, delta, delsq 
/STATES/ q, qdot.qe.qedpt.u, gamma, f ,r,pe 
/FEHD AT/ mass , stif , adamp , bdamp , coxyz , mat , prop , qO , qdotO , 
inertia, scalef 

/IHTEGR/ es,eo,ec,gs , go ,gc,gk, pivot 

/DIKEHS/ ndof ,nact ,nsen ,nesi .mien , jdiag , nnp , nel ,neig 

/COHDAT/ b,hd,hv,fl,f2,ll,12. 


qalpha, qbeta, qalphao, qbetao 
/ELEDOH/ estifm,ndomain,neld,elnum,eldom 
/INTGEK/ f orceid , etype , ix , emat , epr op , pin , id , mask , 

contype , brow , bcol ,hdrow , hdcol , hvros , hvcol , 
b val , hd val , h wal 
/RESEQH/ iadjcy , i count , perm, r Is 


a 


File: acsis.f 
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CsPrograa ACS IS 

C=Purpoee Accelerated CSX Simulation 
C= Author K. Alvin 
C=Date Hay 1990 
CsBlock Fortran 


C — C 

C C 

C prograa ACSXS C 

C C 

C C 

C Purpose: 2nd Order Accelerated CSX Simulation C 

C C 

C C 

C — — C 

prograa ACSXS 

C GET SHARES DATA FILE 


include ’ shared . iac ’ 

C LOCAL VARIABLES 

real«8 t,z(MAXACT) 
integer n, a, runtype, outstrip 
C LOGIC 

call UiPUtCmntype , out skip) 

if (runtype) 100,200,300 

C EXGEHHODE AHALTSXS 

100 continue 

call PREPFEH 
call EIGEHS 
goto 999 

C CSX SIHULATIOH 

200 call PREPFEH 

call PREPCOH 

n a 0 
n a o 

print *, 'Finished Preprocessing . . . starting simulation’ 
print *,’Tiae »*,t0 
CALL ACSISOOT(tO) 

do 250 tatO.tf ,h 

call FQRCES(t*h/2) 

C Predict CSX coupling variables u and gamma 
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if (t . g« . tc) than 
call MEASURE (z) 
call COHTROl(s) 
if (aolag) then 
call SOPELA<I(*) 
call COHTROL(x) 
eadif 
•ad if 

C Structure and observer set ap for parallel execution 

CVS* CHCALL 

do 278 iel,2 

C Integrate Observer Equations 

if CCi .eq. 1) .and. (t .ge. tc)) then 
if (eontype . eq. 0) then 

Call SEC0RDER(aass , stif .adaap.bdamp.f ,go,eo ,qe,qedot 
delta , delsq , j diag , ndof , HAXDCF) 
elseif (eontype .eq. 1) then 
call KFILTER 
endif 

C Integrate Structure Equations 
elseif (i .eq. 2) then 

call SEC0RDER(aass , at if , adaap , bdaatp , f ,go , es , q , qdot , 
delta, delsq, jdiag, ndof , MAID OF) 

endif 

27S continue 

C PRIHT TIME EACH 100 iterations 

a * a + 1 ' 
a « a ♦ 1 

if (n .ge. 100) then 
print », ’Time ■ *,t+h 
n . » 0 
endif 

if (a .go. outskip) then 
call ACSISOUT(t*h) 

erite(24, ’ (40f 12 . 8) * ) t , (z(i) , i»l ,asen) 
a a o 
endif 

2S0 continue 
goto 999 

C TRAHSIEHT RESPOHSE 

300 call PREPFEM 

a a o 
a a o 

print «, ’Finished Preprocessing . . . starting simulation* 
print *,’Tiae a »,t0 
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call ACSISQtJT(tO) 
do 360 tetO,tf,h 
call F0ECES(t*h/2) 
call' ZEHQVECTCgs , ndof) 

call SECQEDE&Cmass , st if , adamp , bdamp , f , gs , e a , q , qdot , 
delta, delsq, jdiag.ndof .HAXDQF) 

C P&XST TIKE EACH 100 iterations 

a ■ a ♦ 1 
b ■ i ♦ 1 

if (a .go. 100) then 
priat e.'Time • ’ ,t*h 
a ■ 0 
eadif 

if (a .go. outskip) tkoa 
call ACSISOOT(teh) 
a ■ 0 
oadif 

3S0 coat inn o 

999 atop 
oad 


File: acsisout.f 


C=Hodule ACSIS0UT 

CaPurpose Vrito doairod output from ACSIS for plotting, otc. 

C=> Author K. Alvin 
C=Date Hay 1990 
C=Block Fortran 

C 

C Subroutine ACSISOUT 

C 

C Purpose: 

C This subroutine outputs formatted displacement and velocity 

C results at a given time for plotting time histories. The 

C desired output variables are defined in outlabel(). 

C 

C 

C Arguments 

C t - time 

C 


subroutine ACSISQUT(t) 

include ' shared. inc* 
real*8 t 
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c LOCAL VARIABLES 

iatagar i 
C LOGIC 

writa(i3, ’ (40212.8) ’ ) t,(q(id(ottlabal(i*10) ,oatlab«l(i))) , 

. i»i,adiaoot) ,(qdot(id(oatlabal(i+30) e oatlabal(i*20)) ) , 

. i»l,aoaloat) 

wrlta<a3,'(40 1*12.8)’) t- > (q a C i d(aat 3a b alCi»lft)«aatlabalCijjJ < 
. i»l „ndiaoat) , (qadot(id(oatlabel(i+30) ,oatlabal(i+20))) , 

. i*l,avalout) 

vrita(2S,’ (40*12.8)’) t,(a(i) ,i=l,nact) ,(gamma(i) ,ial,naea) 

i* (aaiaata) call AHIM0UT(q,id,aap,t,18) 

return 

«od 


File: addstf.f 


— +c 
c 
c 
c 
c 
c 
c 
c 
c 
c 
c 

C 

C ARGUMENTS 

r«al*8 sk(naaq,asaq) , bk(l) 
iatagar ln(18) , jdiag(l) 

C iatagar 1 b( 18)7 jdiag(l), asaq 

C LOCAL ARGUMENTS 

iatagar i, j, k, 1, a 

C ASSEMBLE GLOBAL STIFFNESS AND LOAD ARRAYS 

do 20 j a 1, aaaq 
k « la(j) 

if (k . aq. 0) goto 20 


CsModala ADDSTF 

Caparposa Aaaambla Global ati**naaa matrix 
CaAatbor oho kaooa 
CaUpdata January 1989, by E. Praaoao 
C=Block Fortran 

subroutine ADDSTF(ak,ln,bk, jdiag.aaeq) 

C PURPOSE: 

C THIS SUBROUTINE ASSEMBLES THE ELEMENT STIFFNESS MATRICES 

C INTO THE COMPACTED GLOBAL STIFFNESS VECTOR. 

C 

C ARGUMENTS: 

C ak ELEMENT STIFFNESS MATRIX 

C la LOCATION VECTOR FOR ELEMENT STIFFNESS MATRIX 

C bk COMPACTED GLOBAL STIFFNESS VECTOR 

C jdiag - VECTOR OF DIAGONAL ELEMENT ADDRESSES 

C aaaq - NUMBER OF DEGREES OF FREEDOM PER ELEMENT 
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► 


1 * jdiag(k) - k 

C la jdlag(k*l) - k 

do 10 i » 1, aseq 
n a la(i) 

if(a .gt. k .0E. a .eq. 0) goto 10 

a ■ 1 ♦ a 

bk(a) a bk(a) ♦ ak(i,j) 

10 continue 

20 continue 
C 

return 

aad 

CaEnd Fortran 


File: beam3d.£ 


C=Module BEAM3D 

C=Purpose Conatract 3-d Tiaoahanko baaa aleaant stiffness aad lamped mass 
CaAuthor K. Alvin 
C=Data Nay 1990 
CaBlock Fortran 

sabroatina BEAH3D (n , ni , a j , ak , xyz , eaod , gaod , rho ,araa,ssf2,ssf3 t 
tc jtor,i2,i3, ipia.sk, am) 


C 

C 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 


c 


ARGUMENTS : 

a Elaaant ID Number 

ni Hoda ID Number at End i 

nj Node ID Number at End j 

xyz Node Location Array 

amod Material Elastic Nodalaa (Young's Modulus) 

gmod Material Modulus of Rigidity (Shear Modulus) 

rho Material Mass Density 

area Element Cross-sectional area 

saf2 Shear shape factor in element x2 direction 

ssf3 Shear shape factor in element x3 direction 

jtor Torsional constant J 

12 Area moment of inertia about element x2 axis 

13 Area moment of inertia about element x3 axis 

ipin Pin release codes: 0»Fixed, l=Freed 

(1) Axial 

(2) Torsional 

(3) End A rotation about x2 axis 

(4) End A rotation about x3 axis 

(5) End B rotation about x2 axis 

(6) End B rotation about x3 axis 

sk Element Stiffness Matrix 

sm Element Mass Matrix 

integer n,ni,nj ,nk,ipin(l) 

real*8 xyz(3,l) ,emod t gmod,rho,area,ssf2,ssf3, jtor,i2,i3 
real*8 sk(12,l) ,sm(12,l) 

LOCAL VARIABLES: 

integer i.j 

real*8 dc(3,3) , length , rlength ,kc( 10) ,mc(3) 
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C LOGIC 


C Find El— at Length 

length * O.OdO 
do 10 i * 1,3 

de(i ,i) s xyz(i.nj) - xyz(i,ni) 

Xuagth s length ♦ dc(l,i)**2» 

10 continue 

-length » a qri (length) 
if (length .eq. O.OdO) then 

print *, ’3AH2D: Zero element length; n<* * ,n 

return 

endif 

C Find direction coeinee for xi,x2,x3 element axee 
do IS i*l,3 

dc(l,i) a de( 1,1) /length 
if (nk .eq. 0) then 
de(2,i) » 0.0 
else 

dc(2,i) * xyz(i.nk) - xyzCi.ai) 
endif 

15 continue 

if (nk .eq. 0) dc(2,3) » 1.0 

16 dc(3,l) ■ dc(1.2)*dc(2.3) - dc(2,2)*dc(i,3) 

dc(3,2) ■ dc(2,l)*dc(l,3) - dc(l,l)*dc(2,3) 

dc(3,3) » dc(l »l)*dc(2,2) - dc(2,l)*dc(l,2) 

rlength « eqrt(de(3,l)*»2 ♦ dc(3,2)**2 ♦ dc(3,3)**2) 

if (rlength .ne. 0.) goto 17 

dc(2,2) ■ 1.0 

dc(2,3) ■ 0.0 

goto 16 

17 do 18 i»l,3 

dc(3,i) « dc(3,i) /rlength 

18 continue 

dc(2,l) * dc(3,2)*dc(l,3) - d'c(l,2)*dc(3,3) 

dc(2,2) » dc(l,l)*dc(3,3) - dc(3,l)*dc(l,3) 

dc(2,3) « dc(3,lWc(l,2) - dc(l,l)*dc(3,2) 

C Compute various stiffness constants, accounting for pin codes 

if (ipin(l) . eq. 0) then 
kc(l) ■ are&*emod/length 

else 

kc(l) a O.OdO 
endif 

if (ipin(4) .eq. 0) then 
if (ipin(6) .eq. 0) then 

kc(2) a area*gmod*Bsf 2/length 
kc(6) a i3*emod/length 
ke(7) a kc(2)*length/2.0d0 
kc(9) a kc(7) *length/2 . OdO 
else 

print *,'BEAH3D: Pin code error, x3 direction, el #’,n 
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endii 


alaa 

it (ipia(6) . aq. 0) than 

print * , ’ BEAM3D : Pin coda arror, x3 direction, al #' ,n 

alsa 

kc(2) a O.OdO 
kc(6) * O.OdO 
kc(7) a O.OdO 
kc(9) • O.OdO 
andil 
+ndit 


it (ipin(3) .aq. 0) than 
ii (ipin(S) .aq. 0) than 
kc(3) * are**gBod*aai3/lehgth 
kc(5) a i2*eaod/length 
kc(8) a kc(3)*length/2.0d0 
kc(lO) a kc(8)«length/2.0d0 
alaa 

print *,’BEAH3D: Pin coda arror, x2 direction, al t ’ ,n 
eadii 

alaa 

ii (ipin(5) .aq. 0) than 

print *, ’BEAH3D: Pin coda error, x2 direction, al 8* ,n 

alaa 

kc(3) • O.OdO 
kc(5) a O.OdO 
kc(8) a O.OdO 
kc(lO) a O.OdO 
endif 
endif 

ii (ipin(2) .aq. 0) than 
kc(4) ■ jtor*gBod/length 

alaa 

kc(4) » O.OdO 
endi t 

mc(l) a area«rho*laagth/2 . OdO 
mc(2) a i2*rho*length/2.0d0 
mc(3) a i3*rho*length/2.0d0 


ak(l.l) 

sk(l,2) 

ak(l,3) 

sk(l,4) 

akCl.S) 

ak(l,6) 

ak(l,7) 

ak(l,8) 

ak(l,9) 

ak(l,10) 

ak(l.U) 

ak(l,12) 

sk(2,2) 


» kc(l)«dc(l,l)*dc(l ,1) 
kc(3)*dc(3,l)*dc(3,l) 
» kc(i)*dc(l,l)*dc(l,2) 
kc(3)*dc(3,l)*dc(3,2) 
a kc(l)*dc(l ,l)*dc(l,3) 
kc(3)*dc(3,l)*dc(3,3) 
» kc(7)*dc(2, l)*dc(3,l) 
a kc(7)*dcC2,i)*dc(3,2> 
a kc(7)*dc(2,l)*dc(3,3) 
a -ak(l,l) 
a -ak(l,2) 
a -ak(l,3) 

* ak(l,4) 
a ak(l.S) 
a ak(l,8) 

a kc(l)*dc(l,2)*dc(l,2) 
kc(3)*dc(3,2)*dc(3 ,2) 


♦ kc(2)*dc(2 1 ,l)*dc(2,i) * 

♦ kc(2)*dc(2,l)*dc(2,2) ♦ 

♦ kc(2)*dc(2,l)*dc(2,3) ♦ 

- kc(8Wc(3,l)*dc(2,l) 

- kc(8)*dc(3,i)*dc(2,2) 

- kc(8)«dc(3,l)*dc(2,3) 


♦ kc(2)*dc(2 ,2)»dc(2,2) ♦ 
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ak(2»3) ■ ke(l)*dc(l p2)*dc(l ,3) * kc(2)*dc(2,2)*dc(2.3) ♦ 

kc(3)*dc(3,2)*de(3,3) 

ak(2,4) » kc(7)*dc(2,2)*de(3,l) - kc(8)*dc(3,2)*dc(2,l) 

sk(2,5) « kc(7)*dc(2.2)*de(3,2) - kc(8)*de(3,2)*de(2,2) 

#k(2,8) » ke(7)*dc(2,2)adc(3,3) - kc(8)*de(3.2Wc(2,3) 

#k(2,7) a -sk(i,2) 
ak(2,8) a =sk(2,2) 
ak(2,9) a ~ak<2,3) 

«k(2,10) > ak(2,4) 

ak(2.1i) a ak(2,8) 

»k<3,12) « ak(2,6) 

ak(3,3) ■ kc(i)*dc(l,3)*dc(i,3) <• ke(2)*dc(2,3)*dc(2,3) ♦ 

ke(3)«de(3,3)«de(3,3) 

ak(3,4) a kc(7)*dc(2,3)*dc(3 ,1) - kc(8)*dc(3,3)*dc(2,l) 

sk(3 ,5) a ke(7)*de(2,3)*dc(3,2) - kc(8)*dc(3,3)*dc(2,2) 

ak(3,8) a kc<7)*dc(2,3)*dc(3,3) - kc(8)*dc(3,3)*dc(2,3) 

ak(3,7) » -ak(l ,3) 

sk(3,8) « -ak(2,3) 

ak(3,9) a -sk(3,3) 

sk(3p J.0) a ak(3,4) 

sk(3,ll) a sk(3,S) 

ak(3,12) a S k(3,6) 

8k(4,4) a kc(4)*dc(l ,l)*dc(l ,l)+’(kc(10)+kc(5)) , i'dc(2 ,l)*dc(2 , l) 

♦ (kc(9)*kc(6))*dc(3 , l)*dc(3,l) 

sk(4 ,5) a kc(4)ade(i,i)adca,2)+(kc(10)+kc(S))*dc(2,l)*dc(2,2) 

♦ (kc(9)+kc(6))*dc(3 ,l)*dc(3,2) 

ak(4,6) a kc(4)*dc(l,l)*dc(l,3)*(kc(10)*kc(5))*dc(2.1)*dc(2,3) 

♦ (kc(9)*kc(6))*dc(3,l)*dc(3,3) 

■k(4,7) a -«k(l,4) 

sk(4,8) » -sk(2,4) 

ak(4,9) a -ak(3,4) 

sk(4,10) a -ke(4)*dc(l,l)adc(l,l)*(kc(10)-kcC5))*dc(2,l)*dc(2,l) 

♦ (kc(9)=>kc(6))*dc(3,l)-*dc(3,l) 

sk(4,li) a -kc(4)*dc(l,l)»dc(l,2Wke(10)-kc(5))*dc(2,l)*dc(2,2) 

♦ (kc(9)-kc(8))*dc(3,l)*dc(3,2) 

ak(4,12) a -kc(4)adca,l)*dc(l,3)*(kc(10)-kc(5))*dc(2,l)*dc(2,3) 

+ (kc(9)-kc'(6))«dc(3 , l)*dc(3 ,3) 

sk(5,S) » kc(4)adc(l,2)*dc(l,2)+(ke(10)*kc(S))»dc(2,2)*dc(2,2) 

♦ (kc(9)*kc(6))*dc(3 ,2)*dc(3 ,2) 

sk(5,6) a kc(4)*dc(l,2)*dc(l,3)*(kc(10)*kc(5))*dc(2,2)*dc(2,3) 

♦ (kc(9)*kc(6))*dc(3,2)*dc(3,3) 

sk(5,7) a -ekU.S) 

ak(B,8) a -ak(2,S) 

sk(B ,9) a -ak(3,5) 

ak(S,10) a ak(4,ll) 

akCS.ll) a -kc(4)*dc(l,2)*dc(l ,2)*(kc(10)-kc(5))*dc(2,2)*dc(2,2) 

♦ (ke(9)-kc(6))*dc(3 ,2)*dc(3 ,2) 

ak(5,12) a -kc(4)*dc(l ,2)*dc(l ,3)+(kc(10)-kc(5))*dc(2 ,2)*dc(2,3) 

♦ (kc(9)-kc(6))*dc(2,3)*dc(3 ,3) 

ak(S,6) a kc(4)*dc(l,3)*dc(i,3)*(kc(10)*kc(5))*dc(2,3)*dc(2,3) 

♦ (kc(9)+kc(6))*dc(3,3)*dc(3,3) 

ak(6,7) a -ak(l,6) 

ak(6,8) a -ak(2,6) 

ak(6 ,9) a -sk(3,6) 

ak(6,10) a ak(4,12) 
ak(6,ll) a ak(5,12) 

ak(6,12) a -kc(4)*dc(l,3)*dc(l,3)+(kc(i0)-kc(5))*dc(2,3)*dc(2,3) 

♦ (kc(9)-kc(6))*dc(3 ,3)*dc(3 ,3) 
ak(7,7) a ak(l.l) 
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ak(7.8) a ak(l,2) 
ak(7.9) « ak(l,3) 

•k(7,10) - -ak(l,4) 

ak(7,U) ■ -ak(l,5) 
ak(7 , 12) ■ -ak(l,8) 

ak(8„8) ■ ak(2,2) 

ak(8,9) > ak(2,3) 

ak(8,10) > -ak(2,4) 

sk(8,lt) » -ak(2,6) 
ak(8,12) • -ek(2.6) 

ak(9,9) > ak(3,3) 

ak(9.10) » -ak(3,4) 

ak(9,U) a -ak(3,6) 
sk(9,12) > -ak(3 ,6) 

3k(10,10) « sk(4,4) 

ak(lO,ll) > ak(4,6) 
ak(10,12) - ak(4,6) 
akCll.ll) a ak(S.S) 
ak(ll,12) a ak(5,6) 
ak(l2,12) > »k(8,6) 

do 60 i a 1,12 
do SB j • 1,12 
an(i,j) > o.dO 
55 continue 

SO continue 



C Row-sum rotated mass matrix to re-diagoaalize 


C 

C 

C 

C 

C 

C 


C 

C 

C 

C 

C 

c 


sm(l,l) a mc(l) 

am(2,2) a mc(l) 

sm(3,3) a mc(l) 

am(4,4) a mc(2)*(dc(l,l)*dc(l, l)+dc(2,l)*dc(2,l)) ♦' 

»c(3)*(dc(l ,l)*dc(l,l)+dc(3,l)*dc(3 , 1)) ♦ 
mc(2)*(dc(l ,l)*dc(l ,2)*dc(2,l)*dc(2,2) ) ♦ 

. mc(3)a(dc(l;l)*dc(l,2)-*dc(3,l)*dc(3,2)) ♦ 

mc(2)*(dc(l ,l)*dc(l ,3)+dc(2,l)*dc(2 ,3) ) ♦ 

mc(3)*(dc(l ,l)*dc(l,3)+dc(3,l)*dc(3,3)) 
am(4,S) a mc(2)*(dc(l ,l)*dc(l ,2)+dc(2,l)*dc(2 ,2) ) ♦ 

ac(3)*(dc(l,l)*dc(l,2)+dc(3,l)*dc(3,2)) 
am(4,6) a oc(2)e(dc(l,l)*dc(l,3)+dc(2,l)*dc(2,3)) ♦ 

mc(3)*(dc(l,l)*dc(l ,3)+dc(3,l)*dc(3,3)) 
8a(5,4) a oc(2)*(dc(l,l)*dc(l,2)«>dc(2,l)*dc(2,2)) ♦ 

mc(3)*(dc(l,l)*dc(l.2)t'dc(3,l)*dc(3,2)) 
80(5,5) a oc(2)*(dc(l,2)*dc(1.2)+dc(2,2)*dc(2.2)) ♦ 

oc(3)*(dc(l,2)adc(l,2)*dc(3,2)»dc(3,2)) ♦ 

mc(2)*(dc(l,l)*dc(l ,2)*dc(2,l)*dc(2,2)) ♦ 

mc(3)*(dc(l,l)*dc(l,2)+dc(3,l)*dc(3,2)) ♦ 

oc(2)*(dc(l,2)«dc(l ,3)+dc(2,2)*dc(2,3)) ♦ , 

oc(3)*(dc(l ,2)*dc(l ,3)*dc(3 ,2)*dc(3 ,3)) 
ao(5,6) ■ ac(2)*(dc(l ,2)*dc(l,3)+dc(2,2)*dc(2,3)) ♦ 

oc(3)*(dc(l ,2)*dc(l ,3)*dc(3,2)*dc(3 ,3)) 
ao(6,4) » oc(2)*(dc(l,l)adc(l,3)+dc(2,l)*dc(2,3)) ♦ 

mc(3)*(dc(l , l)*dc(l ,3)+dc(3, l)*dc(3,3)) 
ao(6,5) ■ mc(2)*(dc(l ,2)*dc(l ,3)+dc(2,2)*dc(2,3)) ♦ 

oc(3)*(dc(l,2)*dc(l,3)+dc(3,2)*dc(3,3)) 
ao(8,6) a oc(2)*(dc(l,3)*dc(l,3)+dc(2,3)*dc(2,3)) ♦ 

oc(3)*(dc(l,3)*dc(l,3>dc(3,3)*dc(3,3)) ♦ 
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sa<7,7) 
sa(8,8) 
sa(9,9) 
sh(10,10) 
C sa(iO„ii) 

C s»(10,12) 

C sa(ll,l 0) 

«*(11 9 4D 
C s®(ii,12) 

C sb( 12,10) 

C sa(12,li) 

sa(12,12) 


ac(2)*<de(i,l)*dc(i,3)*dc(2,l)*de(2,3)) + 

ac(3)*<dc<l,l)*dcU,3)+dc(3,l)*dc(3 9 3)) ♦ 

■c(2)*(de(l,2)*dc(l,3)+dc(2,2)«dc(2,3)) * 

ac(3)*(dc(i,2)*de(l,3)+dc(3,2)*de(3,3)) 

ae(i) 

ac(i) 

ae(l) 

sa(4,4) 

s»(4„6) 

sa(4,6) 

ss(5,4) 

.sa£&*S) 

sa(f>,6) 

sa(6,4) 

sb(6,S) 

sm(8,6) 


do 100 i»i,12 
do 200 j»l,i-l 
sk(i, j) a sk(j,i) 
200 continue 

100 continue 


return 

end 

C=End Fortran 


File: forces. f 


CaKodule FORCES 

C=Purpose Calculate applied force rector at given time 
C=Autho r K. Alvin 
C=Date Hay 1990 
C=Block Fortran 

C 

C Subroutine FORCES 
C 

C Purpoae : 

C Returns force from stored function at any given time. 

C The forcing functions are hardwired by the user. The 

C function is selectable at program ezcecution using the 

C forcing function 10, which by convention is the statement 

C label used in branching. 

C 

C ; 

c 


subroutine FORCES (time) 

include • shared . inc ' 
real*8 time, pi 

C LOQXC 

pi > 3. 141S926 
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call ZEHOVECKf.ndof) 

101 if (forcaid . aq. 101) thaa 

if (tima .la. .02) thaa 

f (id(2,15)) ■ 100 . a ( 1 . -coa ( 2 . *pi*t ima/ . 02) ) 
a ad if 

102 alaaif (forcaid . aq. 102) thaa 

if (tima .It. .1) than 
f(id(2.96)) ■ 10. 
alaaif (tima .aq. .1) than 
f (id(2,95)) * 0. 

alaaif ((tima .gt. .1) .aad. (tima .It. .2)) thaa 
f(id(2,96)) « -10. 

alsa 

f(id(2,9S)) a o. 
aadif 

103 alaaif (forcaid .aq. 103) thaa 

if (tima .la. .01) thaa 

f (id(l.lS)) a ioo 
aadif 

104 alaaif (forcaid .aq.104) thaa 

if ((tima .gt. 0) .aad. (tima .It. .17)) thaa 
f(id(3,126)) a io 

alaaif ((tima .gt. .17) .aad. (tima .It. 1.0)) thaa 
f(id(3,12S)) a -io 
aadif 

105 alaaif (forcaid .aq. 105) than 

if (tima .la. .01) thaa 

f(id(2,9)) * I00.*(l.-coa(2.*pi*tima/.01)) 
alaaif (tima .la. .02) thaa 

f(id(2,9)) ■ 100 . *( cos (2 . *pi*tiaa/ . 01 )-l . ) 
aadif 

aadif 

do 10 i a i, adof 

f(i) * acalaf * f(i) 

10 coatinaa 

ratarn 

aad 


File: input. f 


C=Hodula IS POT 

C=Puxpoaa Input data paramatara for ACSIS 
C= Author K. Alvin 
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CsDate Hay 1990 
C=Block Fortran 


C 

C 

C 

C 

C 

C 

C 

C 

C 

C 

c 

c 

c 

c 


Subroutine INPUT 


Argument definitions 


runtype 

savin 

runfile 

COBBSnt 

outskip 


10 of analysis run type 

variable to control creation of input file 
variable indicates if run is froa input file 
dummy naae for comment input lines 
nuaber of steps to skip before sending output 


subroutine INPUT (runtype, out skip) 


include ' shared. inc' 
integer runtype, outskip 
character 1 *! savin, ruaf ile, temp 
char act er *48 comment , inpf ile 


C PRINT AND READ START-UP 


print *, ’2nd Order Accelerated CSI Simulation (ACSIS)' 
print * 

print *, 'Please input analysis type:* 
print * 

print «, * 1. Eigenmode Analysis' 

print *, ’ 2. CSX Simulation' 

print *, ' 3. Transient Response' 

print • 

read ■*, runtype 


C RUN OPTIONS AND INPUT FILE SETUP 


runfile ■ 'n' 
if (runtype .It. 0) then 
runtype -1 * runtype 
runfile » 'y' 
endif 

print *, 'Do you oish to save an input file? (y or n)’ 
read 21, savin 

20 format (a32) 

21 format (ai) 

if (savin .eq. * y ’ ) then 

print *, ’Name of save input file? (filename)’ 
read 20, inpf ile 
open( 16 , f ileainpf ile ) 
runtype a -1 *• runtype 
vrite(16 , ’ (i2) 1 ) runtype 
. runtype » -1 * runtype 
vrite(16, ’(al)*) ’n' 

write(16,’(a47)’) '* ACSIS input file,tvo lines above are' 
vrite(16,’(a48)') '* analysis type and save input file. Do' 
vrite(16, ’ (a48) ’ ) ’* not change them by editing this file.’ 
endif 
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ruatype a runtype - 2 

if (runfile .aq. ’y’) than 
do 30 i a 1,4 

road 20, comment 
30 eontinno 

endif 

print *, ’Finite Element Xodal Input File Name (filename)’ 
road 20, femfile 
open (11 ,f ile=f oaf ilo) 
if (savin .aq. ’y’) than 

writo(18, ’ (a47) ’ ) ’* Finite element input file?(f ilename) * 
writoClfl, ’ (a32) ’ ) f oof ilo 
endif , 

if (runtypo) 100,200,300 

C EIGENHODE INPUTS 


100 print a, ’Busbar of modos desired: ’ 

if (runfile .aq. 'y') road 20, comaant 
road *, neig 

if (runfila .aq. *y') read 20, comment 
print a, 'Output File Hama:’ 
read 20, outfile 
open(13,file a outfile) 
if (savin .aq. ’y') than 

vrita(16, ’(a36)’) ’a Humber of modes desired? 1 
vrite(l6,’(i4)’) naig 

vrite(16,’(a33)’) ’a Output file?(f ilename) ’ 
writ a (16 , ’ (&32) ’ ) ontfila 
andif 


call READ FEB 


goto 1000 
C CSI INPUTS 

200 print a ,» Controller Definition File Name:' 
if (runfile .aq. ’y') read 20, comment 
read 20, confile 
open(l2 ,file s conf ile) 

print a, ’Please input type of control:’ 
print * 

print *, ’ 1. Full State Feedback' 

print a, » 2. Luenbarger Observer (1.1=0) ’ 

print a, » 3. Kalman Filter' 

print * 

if (runfila .aq. *y’) read 20, comment 
read *, contype 
contype a contype - 2 

print *, ’Initial time, final time, control-on time, step size: 
if (runfile .eq. ’y’) read 20, comment 
read a, tO,tf,tc,h 

print a, 'Forcing function ID, scale factor, damping coeff- a.b 

if (runfila .aq. *y*) read 20, comment 

read «, forceid.scalef ,adamp,bdamp 

print *, ’Phase lag fix?(y or n):' 

if (runfile .aq. ’y’) read 20, comment 
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read 21, tamp 

if (tastp ,eq. 'y' ) nolag * .true, 
if (temp . eq. 's') nolag ■ .false, 
print a, 'Gaia scale factors (4 total):* 
if (ruafile .eq. *y’) read 20, cement 
read * , qalph* 0 qbeta , qalphao , qbetao 
if (savin .eq. *y’) then 

write(16, *(a42)’) ** Costroller file name?(f ilename) * 
write(16,’(a32)*) eosfils 

write(16,'(a42)') '• Please iapat type of costrol: * 

write(16, *(iiO)') costype * 2 

write(16,*(a46)’) ’* Xaitial.final.eontrol-oa.step size?* 
write( 16 , » (4f 14 . 8) * ) t0,tf,tc,h 

vritq(16,’(a49)*) ** Forcing fsactioa, scale f, damping a.b? 
write(16,*(i4,flS.6,2fl2.8) *) forceid,scalef ,adamp,bdamp 
write(16, ' (a32) * ) •* Phase lag fix?(y or a)’ 
if (aolag) write(16,'(al)*) *y* 
if (.aot. aolag) write(16, *(al)*) 'a' 
write(16,*(a40)') *« Gaia seals factors (4 total)?* 
writ e( 16 , ’ (4f 14 . 8) * ) qalpha , qbeta , qalphao , qbetao 
eadif 

call READFEH 
goto 999 

C TRANSIENT RESPONSE INPUTS 


300 print *, 'Initial time, final time, step size:' 
if (runfile .eq. *y*) read 20, comeat 
read *, tO.tf.b 

print *, 'Forcing fsactioa ID, scale factor, damping coaff- a.b 
if (raafile .eq. *y') read 20, comeat 
read «, f orceid , scalef ,adamp,bdamp 
if (savia .eq. *y*) then 

write(16, *(a48)') ’* Initial, final, step size?* 
write( 16 , * (3f 14 . 8) ’ ) tO.tf.h 

write(16,*(a49)*) ’* Forcing function, scale f, damping a.b? 
write(16, *(i4,flS.6,2fl2.8) *) f orceid, scalef ,adamp,bdamp 
eadif 

call REAOFEM 


goto 999 

C OUTPUT OPTIONS 

999 print *, 'Output File Name:* 

if (ruafile .eq. *y*) read 20, comeat 
read 20, ostfile 
open( 13 , f ile»outf ile ) 

print *, ’Number of displacement results to output (max 10):' 
if (ruafile .eq. *y') read 20, comeat 
read *,ndisout 
do 600 i s l,ndisout 

print * , ’ Input node *, dot tor displacement output#' ,i 
read *,outlabel(i),outlabel(i+10) 

500 continue 

print *, 'Number of velocity results to output (mar 10):' 
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if (runfila .aq. * jr' ) road 20, comment 
read * ,nvelout 
do 600 iai.nvelout 

print * , ' Input nod* #, dof for velocity outputs', i 
road * ,outlabel(i+20) ,outlabel(i+30) 

600 eontinuo 

print *,'Send output ovory boo aany stops?' 
if (ruafilo .aq. *y') road 20, comment 
road * , outskip 
if (savia .oq. ’y‘) tboa 

vrito(16, '(a38)') '* Output fils aaao?(filoaaiBo) ' 
write(16,'(a32)*) outfilo 

writs (16, ’(a42) *) '» Humber of displacement outputs?' 
vrite(16, *(i4)*) ndisout 
do 660 i*l .adisout 

write(16,'(2i8)') outlabal(i),outlab*l(i*10) 

660 eontinuo 

srito(16, ’(a38)') '* Huabor of volocity outputs?' 
mrite(16, *0L4)O avolout 
do 660 ial,nvelout 

writo(16, * (2i8) ' ) outlab*l(i*20) , outlabal(i*30) 

660 continue 

orito(16, *(a44)*) ’* Send output ovary bow many stepj?’ 
srite(16, ’ (i3) ’) outskip 
ondif 

C AHIHATIOH 0PTI0H 

print * , * Animation Output? (y or n) : ’ 
if (runfila .aq. 'y') road 20, commant 
road 21, tamp 

if (tamp .aq. *y') animats « .true, 
if (tamp .oq. 'n') animats » .falsa, 
if (animats) tban 

print*, 'Animation fils name (filename) ' 
if (runfila .aq. 7 'y') road 20, comment 
road 20, animfilo 
open (15, file*animfilo) 
andif 

if (savin .aq. *y') tboa 

vrite(16,*(a41)*) ’• Sand animation output?(y or n)' 
if (animats) vrito(16,'(al)') 'y' 
if (.not. animats) vrito(l6,’(al)') 'n' 
if (animats) tboa 

vrito(16, ' (a31) ') ’* Animation file name?' 
vrito(16,'(a32)*) animfilo 
andif 
andif 

delta * b/2. 
dolsq a dolta**2 

1000 return 
end 


File: pmvmul.f 


61 





CsModale PMVHUL 
CsAntkor g. Airis, 
CsDate May 1990 
C»Bleek Fortran 


C 

C Subroutine PKVHUL 
C 

C Purpose: 

C T&ie subroutine multiplies a matrix ia vector lorn 

0 aad a vector. 

C 
C 
C 
C 
C 
C 
e 
c 
c 
c 


Arguments 

a - matrix ia vector form 

b - vactor 

c - reaslt vactor 

aeq - order of vactor aad square matrix 

jdiag - array of diagonal addresses for a 


C aabroatiaa PHVHDLCa, jdiag, b.neq.c) 


recursive aabroatiaa PHVH0L(a, jdiag, b, aaq, c) 

raal*8 a(l), b(l), c(l) 
integer jdiag(l), aaq 

do 100 i®l ,aaq 

c(i) » a(jdiag(i))*b(i) 

100 coatiaaa 


do 200 ia2,aaq 

do 300 jajdiag(i-l)+l,jdiag(i)-l 
k » jdiag(i) - j 
c(i) ■ c(i) ♦ a(j)*b(i-k) 

300 coatiaaa 

200 coatiaaa 

do 250 ia2,aaq 

do 400 jajdiag(i**l)+l , jdiag(i)-l 
k ■ jdiag(i) - j 
c(i-k) ■ c(i-k) ♦ a(j)ab(i) 

400 coatiaaa 

250 coatiaaa 

ratora 

aad 

C 

C 

C 

C Subroutiaa PMVHAD 

C 

C Porpoaa : 

C Multiply a. matrix ia vactor form aad a vactor and add the 

C. raaultaat vactor multiplied by a constant to a second vactor 
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multiplied by a second vector 


} 


C 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 


•100 


300 

200 


400 

250 


Arguments 

a 

b - 
e 

factl - 
fact 2 - 

jdiag - 
neq 


matrix ia vactor form 

vactor to ba multiplied with matrix 

rasaltaat aad vactor to ba added 

constant multiplier of matrix and first vactor 

constant multiplier of sacond vactor 

array of diagonal addresses for matrix 

order of vectors aad matrix 


subroutine PMVMAD(a, jdiag, b.neq, factl ,c,fact2) 

recursive subroutine PMVHAD(a, jdiag, b, neq, factl, c,fact2) 

xeal«8 «&)„ b£l},c&), fasti ,fact2 
integer jdiag(l), neq 

do 100 i«l,neq 

c(i) ■ fact2«c(i) ♦ factl*a(jdiag(i))*b(i) 
continue 

do 200 1*2, neq 

do 300 j»jdiag(i-l)*l, jdiag(i)-l 
k - jdiag(i) - j 
c(i) a c(i) ♦ factl*a(j)*b(i-k) 
continue 
continue 

do 250 ia2,neq 

do 400 j»jdiag(i-l)+l,jdiag<i)-l 
k a jdiag(i) - j 

c(i-k) » c(i-k) ♦ factl*a( j)*b(i) 
continue 
continue 

return 

end 


\ 


File: prepfem.f 


CaHodule PREPFEH 

C=Purpose Preprocess Structure Finite Element module for ACSIS 
C= Author K. Alvin 
CaDate May 1990 
C=Block Fortran 

C 

C 

C Subroutine PREPFEH 
C 

C Purpose: 

C This subroutine prepares the finite element mass, 

C stiffness, and S matrices in reduced profile vector form 
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c 

c 

c 

c 

c 

c 

c 

e 

c 

c 

c 


Local variable*: 

sk Elaaaat Stiffness aatrir 

an Element Xui Matrix 

3a Local/Global DOF Mapping vector 

nseq Humber of olaaeat degrees of freedom 

ea,ep Material aad Propety id * for eleaeat 

subroutine PREPFE* 


include ’ shared. inc' 


C LOCAL VARIABLES 

paraaetor (MAXSEQ*24) 

ml *8 akOUlSEQ.MAZSEgj.saauZSSq.JUJSSQ^^e^e 
integer la(MAXSEQ) ,nseq,em,ep 

call REBUM 

C Set up akyliaa storage profile for global matrices 

call PR0FXLE(ix,id,jdiag,nnp,ael,4,6,mlen,ndof .mask) 

C Perform automatic domain decomposition 
call DOMDEC 

C Check six* of skyline profile against storage limitation 
if (alen .gt. MAXMLEB) then 

print*, 'PREPFEM: error, global matrix exceeded max. size' 
endif 

C Zero Global Matrices prior to assembly 

call ZEROVECT(stif ,mlen) 
call ZEROVECT(maas .mien) 

C ASSEMBLE EACH ELEMEBT MASS AMD STIFFBESS 

do 100 n*l,nel 

do 20 k*l,4 
j*ix(k,n) 

if ((etype(n).eq.l).and.(k.gt.2)) j » 0 
do 30 i*l,6 
kk=6*(k-l) ♦ i 
if (j .ne. 0) then 
la(kk) * id(i,j) 
else 

la(kk) * 0 
endif 

30 continue 

20 continue 
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if (atypa(a) . aq. 1) tkaa 
aaaq ■ 12 
•a ■ aaat(a) 

«p * aprop(a) 

call BEAM3D(n,ix(l,a) ,ix(2,a) t ix(3,n} ,coxyx,mat(l,8m) , 
aat(2, an) ,mat(3,aa) , prop(i, ap) ,prop(S,ep) ,prop(6,ep) , 

. prop (2 , ap) ,prop(3, ap) ,prop(4 , ap) , pia( 1 ,a) , sk , an) 

alaaif (ix(l,a) .a«. 0) tkaa 

print*, 'PREPFEH.-Eleaaat t jrpa not found, a 3 ' , a, ^types’ ,atype(a) 
aadif 

C ADD ELEMENT TO GLOBAL MASS AMD STIFFNESS 

call ADDSTF(ak,lm,atif , jdiag.aaaq) 
call ADDSTF(an,la,naaa, jdiag.aaaq) 

C SAVE THE ELEHEHT STIFFNESS FOE E-8Y-E COMPUTATIONS 

call SAYESE(ak*a,JU«q) 

100 coatiaaa 

C ADD LUMPED INERTIAS TO GLOBAL MASS? 

do 12S i=l,nnp s 

do 130 j«l,fl 

it (id(j,i) -«q- 0) goto 130 
k»jdiag(id(j,i)) 

maas(k) > aaaa(k) ♦ iaartia(j.i) 

130 coatiaaa 

12S coatiaaa 

C ASSEMBLE AND FACTORIZE aa (S MATRIX) 

me a 1. ♦ dalta«adamp 
kc » dalta*bdamp ♦ dalaq 
do 200 i 3 l,nlaa 

aa(i) 3 mc*maaa(i) ♦ kc*atif(i) 

200 coatiaaa 

♦ 

call SQLVER(aa,gn,jdiag,adof ,1) 

C INITIALIZE DISPLACEMENT AND VELOCITY VECTORS 

do 300 i 3 1 , aap 
do 3S0 j 3 1,6 

if (id(j,i) -aa. 0) thaa 
q(id( j ,i)) 3 qO(j,i) 
qdot(id( j ,i)) 3 qdotO(j.i) 
aadif 

3S0 coatiaaa 

300 coatiaaa 

ratora 

aad 

aabroatina SAVESX(ak,a,aaaq) 
iaclada ’ akarad.iac’ 
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real*8 sk(aseq,i) 
integer a.aseq 

k*0 

do 10 j»i„aseq 
do 20 1*1, j 
k»kel 

*eti£®OE,a)*skd,j) 

20 eoatiaae 

10 co&tiaae 

return 

•ad 

subroutine DOHDEC 

iaclado ’ shared. iac' 

logical aehk,adehk(HAXHODE,)UXDOK) 
integer ados 

do 10 j-l.lttXDOH 
a«ld(j)«0 
do 20 i*i,aap 

adebkd , j)» .false . 

20 eoatiaae 

10 eoatiaae 

ados&iaeO 

do 100 aal.nel 

ndo®*0 

nchJc»0 

do 200 vhile (achk.eq.O) 
ndom»ndom+l 

if (adoa.gt.adoaaia)- adoaaiaaadon 
aehk*l 

if (adchk(iz(l,a) .adorn/) nchk=0 
if (adcltic(ix(2,a),adom)) achk*0 
if (achk.eq.l) then 
•ldo»(n)*adom 
adcUcdzd.a) , adorn)*. true. 
adehk(ix(2,a) ,adon)».tnie. 
end if 

200 eoatiaae 

ueld(adon) *aeld (adorn) ♦! 
elaan(aeld(adoa) , adorn) *a 

100 eoatiaae 

retara 

ead 


File: profile.f 
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•» * 


C=Module PROFILE 

CsPurpoee Compute the number o f equations and set profile for K 
C»Autkor Bob Taylor 
C=Date wbo knove 

CsOpdate January 1939 by E. Pramono 
C= Block Fortran 


subroutine PROFILE ( ix,id, jdiag, nnp , nel , uen , adof , and , neq , mask) 

C+ ♦C 

C PURPOSE: C 

C THIS SUBROUT IHE COMPUTES THE HUMBER OF EQUATIONS REQUIRED C 

C TO SOLVE THE PROBLEM BT EL IH HUT INS RESTRAINED DEGREES OF C 

C FREEDOM FROM THE SYSTEM OF EQUATIOHS. KNOWING THE EQUATION C 

C NUMBERS COORESPOND IHG TO THE HODAL DEGREES OF FREEDOM, THE C 

C DIAGONAL ELEMENT LOCATIONS CAN BE COMPUTED FOR STORING THE C 

C GLOBAL STIFFNES MATRIX IH COMPACTED VECTOR FORM. C 


C 

C ARGUMENTS 
C 

integer ix(aea,l), id(ndof,l), jdiag(l) 
integer nnp, nel, nad, aeq, mask(l) 

C integer nnp, nel, nen, adof, nad, aeq, mask(l) 
C 

C .LOCAL ARGUMENTS 
C 

integer i, j. k. 1, a, n, jl, kl, 11, ml 
C 

C SET UP EQUATION NUMBERS 
C 

aeq a 0 

do 30 n > 1, nnp 
do 20 a a 1, adof 
j » id(a,mask(n)) 
if (j . eq. 1) goto 10 
neq » aeq + 1 
id(a,mask(n)) ■ neq 
jdiag(neq) a 0 
goto 20 

10 id(a,maak(n)) a o 

20 continue 

30 continue 
C 
C 

C COMPUTE COLUMN HEIGHTS 
C 

do 80 n a i, ael 
do 70 a ■ 1, nea 
ml a ix(m,a) 
if (ml .le. 0) goto 70 
do 60 1 ■ 1, adof 
11 a id(l,al) 
if (11 .eq. 0) goto 60 
do SO k ■ a, nea 
kl > ix(k,a) 
if (kl .le. 0) goto SO 
do 40 j si, adof 
jl » id(j.kl) 
if (jl .eq. 0) goto 40 
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i a HAXO(li,jl) 

jdiag(i) = MAXO( jdiagCi) , XABS(H“ji)) 

40 continue 

SO continue 

SO continue 

70 continue 

80 eeatiaue 

C 
C 

C COMPOTE DXAGOHAL PO UTTERS 
C 

aad * I 
jdiag(i) ■ 1 
if (neq . cq. 1) return 
do 90 i * 2, neq 

jdiag(n) » jdiag(a) ♦ jdiag(a-l) ♦ i 
90 continue 

aad a jdiag(neq) 

C 

return 

•ad 

C=Ead Fortran 


File: read.f 


CaRodule READ 
C= Author K. Alvin 
C^Date May 1990 
CaBlock Fortran 

C — 

c 

C Subroutine READFEM 

C 

C Purpoa* : 

C This aubroutin* roads fh« data file for the finite 

C element model. 

C 

C 

C Arguments 

C ctype - stores code for type of line 

C 

subroutine READFEM 
C GLOBALS 

include 'shared.iac' 

C LOCALS 

integer j ,a,ctype,GETTYPE 
chaxacter*132 aline 
real*8 in 

C INITIALIZE SIZE OF PROBLEM 
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sap a 0 
aal * 0 
ndof a 0 
adoaaia a 0 


C IDENTIFY CARD TYPE ABO ASSIGN DJPUT 


10 

100 

150 


raad(U,1000,aada9999) alias 
ctypa » GETTYPE( alias) 
if (ctypa) 10,10,150 


if (aliaa(l:4) .sq. 
if (aliaa(l:4) .sq. 
if (alias(l:4) .sq. 
if (alias(l:4) .sq. 
if (alias(l:4) .sq. 
if (alias(l:4) .sq. 
if (aliaa(l:4) .sq. 
if (alias(l:4) .sq. 
if (alia*(l-:4) .sq. 
if (alias(l:4) .sq. 


* HODS') goto 200 
’TOPO’) goto 300 
'ATTR') goto 400 
'MATE’) goto 500 
* PROP ’ ) goto 800 
’Fill’) goto 700 
'nrrr) goto soo 

•INER’) goto 900 . 
>EH9 >) goto 10 
’MESH’) goto 10 


priat a.'REAOFEH: Uaracogaizad cud typo; 
goto 10 


’ ,aline(l:4) 


C READ MODES 


200 road(ll,1000,sada9999) alias 

ctypa a GETTYPE(alias) 
if (ctypa) 200,250,100 

250 rsad(alias ,*) a, (coxyz(j ,a) , j»l ,3) 
if (a .gt. aap) nap a a 
goto 200 

C READ TOPOLOGY 

300 r aad( 11 , 1000 , aad»9999 ) alias 
ctypa » GETTYPE(aliaa) 
if (ctypa) 300,350,100 . 

350 rsad(aliaa,«) a.atypa(a) ,(ix(j ,a) , j»l ,4) 
if (a .gt. asl) aal a a 
goto 300 

C READ ATTRIBUTES 

400 raad(ll,1000,aada9999) alias 

ctypa ■ GETTYPE(aliaa) 
if (ctypa) 400,460,100 

450 raad(alias,«) a,aoat(a) , eprop(a) , (pia(j ,a) , j=l ,6) 
if (aldon(a).gt.adooaia) adoaaia s eldoa(a) 
goto 400 

C READ MATERIAL 

500 raad(ll,1000,aada9999) aliaa 

ctypa a GETTYPE(aliaa) 
if (ctypa) 600,550,100 

550 raad(aliaa,«) a,(n>at(j,a),jai,3) 
goto 500 

C READ PROPERTIES 
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600 r®ad(ll,1000,«nds9999) alia® 

ctyp® » GETTYPE( aline) 
if (etype) 600,660,100 
660 r®ad(aliac,*) a,(prop(j,a),j»l,6) 
g ota 600 

C BEAD Film 

700 r®ad(ll,1000,®ndo9999) alia* 

etype » OETTYPE(allae) 
if (etype) 700,760,100 
760 r®ad(alia«,*) n,(id(j ,u) , j«l ,6) 
goto 700 

C BEAD ramiL C0HDZTI0HS 

800 read(ll,1000,«nd<>9999) alia® 

etype • SEmPE<a2ia«) 
if (ctyp®) 800,860,100 
860 r«ad(alin®,*) n,j,q0(j,n),qdot0(j,a) 
goto 800 

C READ IBERTIA 

900 read(ll,1000,endo9999) alia® 

etype a GETTTPE(alia®) 
if (ctyp®) 900,960,100 
950 read(alia«,*) a.j.ia 

in«rtia( j ,n)»inertia( j ,a)*in 
goto 900 

1000 foraat(al32) 

9999 coatiau® 

return 

and 

C 

C 

C Snbroatia® READC0H 

C 

C Parpoa®: 

C This anbrontia® reads th® actuator and sensor 

C locations and th® gains for th® control system 

C 

C 

C Arguments 


c 

c 

ctyp® 

* 

stores cod® for type of line 

c 

c 

Abbreviations 

c 

HACT 

- 

nnmbor of actuators 

c 

HSEH 

- 

number of sensors 

c 

BKAT 

- 

locations of actuators 

c 

HD HA 

- 

array of displacement sensor locations 

c 

HVHA 

- 

array of velocity sensor locations 

c 

F1GA 

— 

control gain matrix 
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C F2GA ~ control gain matrix 

C L1GA - atata aatimator filter gain matrix 

C L2GA - atata aatimator Iiltar gain matrix 

anbrontina READCOH 

induda ’ ahared . inc ' 

C LOCALS 

raal*8 val 1 

intagar j ,n, ctypa, GETTYPE 
character* 132 alina , 

bval a 0 
hdval a 0 
hvval » 0 
hdbv&l » 0 
hvbval ■ 0 

C IDENTIFY CARO TYPE AND ASSIGN INPOT 

10 raad(12,1000,end*9999) alina 

100 ctypa » GETTYPE(alina) 

if (ctypa) 10,10,150 

150 if (aline(l : 4) .aq. 'NACT') goto 200 

if ( alina(l : 4) .aq. 1 NSEN 1 ) goto 300 
if (aline(l:4) .aq. 'BRAT') goto 400 
if (aline(l : 4) .aq. ’NONA*) goto 500 
if (alina(l:4) .aq. ’HVRA* ) goto 600 
if (alina(l:4) .aq. 'F1GA') goto 700 
if (alina(l:4) .aq. >F2GA>) goto 800 
if (alina(l:4) .aq. *L2GA*) goto 900 
if (alina(l:4) .aq. 'LiGA') goto 1100 
if (alina(l:4) .aq. 'END ') goto 10 

print «,'READC0N: Unrecognized card type; ' ,alina(l:4) 
goto 10 

C READ INPUT CARDS 

200 raad(12,1000,end;>9999) aline 

ctypa a GETTYPEC alina) 
if (ctypa) 200,250,100 
250 r«ad(aline,*) nact 
goto 200 

300 raad(12,1000,and a 9999) alina 

ctypa a GETTYPEC alina) 
if (ctypa) 300,350,100 
350 readCaline,*) naan 
goto 300 

400 raad(12,1000,and39999) alina 

ctypa a GETTYPE(alina) 
if (ctypa) 400,450,100 
460 read(aline,«) i,j,n,val 
bval 3 bval * 1 
b(bval) 3 vai 
brov(bval) a id(j,i) 
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bcol(bval) a a 
goto 400 

600 raad(12,1000,aad«9999) alia* 

ctypa ■ GETTTPE(aliaa) 

. if (ctypa) 500,660,100 
660 raad(alias,«) i,J,a,»al 
bdval a bdval ♦ 1 
bd(bdval) a v&l 
hd7o«(bdTal) a a 
bdeol(bdval) ■ ld(j,i) 
goto 600 

600 road(12 , 1000 , sad»9999) alia* 
ctypa ■ GETTYPEC alias) 
if (etypa) 600,680,100 
650 rsad(aliaa,«) i,j,a,val 
bwal ■ hwal ♦ 1 
hv(bvval) a -ral 
hvroo(hwal) » a 
bveol(brwal) a id(j,i) 
goto 600 

700 rsad(12,1000,sada9999) alia* 

ctypa a GETTYPEC alias) 
if (ctypa) 700,750,100 
760 raad(alias,*) i,j,a,val 

fl(a,id(j,i)) » qalpha*val 
goto 700 

800 rsad(12, 1000 . snd=9999) alias 
ctypa a GETTYPE(aliaa) 
if (ctypa) 800,860,100 
860 raad(aliaa,«) i,j,a,val 

f2(a,id(j ,i)) a qbata*val 
goto 800 

900 raad( 12 , 1000 , end a 9999 ) alias 
ctypa * GETTYPEC alias) 
if (ctypa) 900,960,100 
950 rsad(aliaa,«) i,j,a,val 

12(id(j ,i) ,a) a qbatao*val 
goto 900 

1100 raad(12,1000,sada9999) alias 

ctypa a GETTYPE(aliaa) 
if (ctypa) 1100,1150,100 
1150 raad(aliaa,«) i,j,n,val 

llCidCj ,i) ,a) a qalpbao«val 
goto 1100 

1000 fonaat(al32) 

9999 coatiaaa 
ratnra 
sad 

C 

C 
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C Function GETTYPE 

C 

C Purpose : 

C This function identifies whether a line is a character 

C input, data input or comment . 

C 

C 

C 

function GETTYPE(string) 

C GLOBALS 

character*132 string 
C LOCALS 

integer GETTYPE, ctype( 10) 
character*! head(lO) 

data head /•!*,*•», 't* ,'$* ,»X» , **»,»•* ,»C» ,'c* .» */ 
data etype ,-1,-1 ,-1,-1, -1,-1 ,0/ 

C LOGIC 

GETTYPE* 1 
do 100 i*l,10 

if (string(l:l) .eq. head(i)) GETTYPE=ctype(i) 

100 continue 

return 

end 


File: solver. f 


C*Kodule SOLVER 

C=Purpoae Solves the system of linear symmetric equations 

C= Author who knows 

C=Date 


C=0pdata January 1989 by E. Pramono 
C=Block Fortran 

C SUBROUTINE SOLVER ( BK , BR , JD I AG , HEQ , IFLAG ) 

recursive SUBROUTINE SOLVER ( BK , BR , JD I AG , NEQ , IFLAG ) 


C+ 

C 

C 

C 

C 

C 

C 

C 

C 

C 

C 

C 

C 

C 


PURPOSE: 

THIS SUBROUTINE SOLVES THE SYSTEM OF LINEAR SYMMETRIC 
EQUATIONS IN VECTOR FORM USING THE CROUT REDUCTION 
METHOD. 

ARGUMENTS: 

BK - GLOBAL STIFFNESS EQUATIONS IN VECTOR FORM 

BR - GLOBAL LOAD VECTOR 

JDIAG - LOCATION VECTOR FOR DIAGONALS IN [BK] 

NEQ - NUMBER OF EQUATIONS 

IFLAG - FLAG INDICATING WHICH FUNCTION IS TO BE PERFORMED 

1 -> FORWARD REDUCTION 

2 -> BACKWARD SUBSTITUTION 


+C 

C 

C 

C 

C 

C 

C 

C 

C 

C 

C 

c 


C 

c 
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Q ARGUMENTS 


BI<1), BE(1) 

INTEGER JDIAG(l), NEQ, I FLAG 
C LOCAL VARIABLES 

REAL*8 ZERO, EZERO, TOL, DAVAL, DOT, D, HDD, DO 

integer ldflag, jr, j, jd. jh, is, ie, k, jdt 

INTEGER JJ, ID, I, IB, IB 


C 

C 

C 

C 


c 

c 

c 

c 


c 

10 


c 

c 

c 

c 


c 

c 

c 

c 


20 

c 

c 

c 


JJ a 6 


NEW PARAMETERS 


2ERS • 0.000 
EZERO > 0 . 3D- 14 
TOL * O.SD-7 
LDFLAG » 0 


FACTOR BK TO UT*D*U OR REDUCE R 


JR a 0 

do 70 j ■ l, req 
JD a JDIAG(J) 

JH « JD - JR 

IS a J - JH ♦ 2 

IF (JH - 2) 60, 30, 10 

IF (IFLAG .HE. 1) GOTO GO 
IE a J - 1 
K a JR ♦ fi 
ID » JDIAG(IS-l) 


IF DIAGONAL IS ZERO COMPUTE A HORM FOR SINGULARITY TEST 


JDT ■ JDIAG(IE) ♦ 1 

IF (BK(JD) .EO. ZERO .AND. IFLAG .EQ. 1) THEN 
CALL DATEST (BK(JDT), JH-2, DAVAL) 

END IF 


REDUCE ALL EQUATIONS EXCEPT FIRST ROW AND DIAGONAL 


DO 20 I a IS, IE 
IR « ID 
ID a JDIAG(I) 

IH « MIN0(ID-IR-1,I-IS*1) 

IF (IH .GT. 0) BK(K) » BK(K) - DOT(BK(K-IH) , BK(ID-IH), IH) 

K ■ K * 1 
CONTINUE 


REDUCE FIRST ROW AND DIAGONAL 
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o o o o noon noon 


C 

30 


40 


IF (IFLAG .HE. 1) GOTO 50 
IE ■ JE ♦ 1 
IE » JO - 1 
K ■ J - JO 
DO a BK(JD) 

DO 40 I ■ IE, IE 
ID ■ JDIAG(K«T) 

IF (BKCID) .EQ. 0.0) GOTO 40 
0 » BK(I) 

BK(I) « BK(I)/BK(ID) 

BK(JD) ■ BK(JD) - D*BK(I) 
CONTINUE 


CHECK FOE POSSIBLE EBEOES AND PRINT WARNINGS 


ROD a BK(JD) 

IF (DABS (EDO) -LX. XOL*DABS(DB)) 

IF (DO .LT. ZERO .AND. ROD .GT. ZERO) 
IF (DO .GT. ZERO .AHD . ROD .LT. ZERO) 
IF (DABS(RDD) .LT. EZERQ) 


WRITE (JJ.2000) J 
WRITE (JJ ,2001) J 
WRITE (JJ ,2001) J 
WRITE (JJ ,2002) J 


COMPLETE RANK TEST FOR A ZERO DIAGONAL TEST 


IF (DD .EQ. ZERO .AHD. JH .GT. 0) THEM 

IF (DABS(RDD) .LT. TOL*DAVAL) WRITE (JJ.2003) J 

END IF 


REDUCE RIGHT HAND SIDE 

50 IF (IFLAG .Eq. 2) BR(J) » BR(J) - DOT(BK(JR*l) , BR(IS-l), JH-1) 

60 JR a JD 

70 CONTINUE 

IF (IFLAG .HE. 2) RETURN. 

C 

C DIVIDE BY DIAGONAL TERMS 

00 80 I a i ( HEq 
ID a JDIAG(I) 

IF (BK(ID) .HE. 0.0) BR(I) a BR(I)/BK(ID) 

IF (BR(I) .HE. ZERO) LDFLAG a 1 
80 CONTINUE 

C 

C 

C CHECK FOE ZERO LOAD VECTOR 

IF (LDFLAG .EQ. 0) WRITE(JJ,2004) 

C 

C BACK SUBSTITUTE 4 

J a HEQ 
JD a JDIAG(J) 

90 Da BR(J) 

J a J - 1 
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IF (J .LE. 0) RETURN 
JR > JDXAG(J) 

IF (JO - JR .LE. 1) GOTO 110 
XS»J-JD*JR*2 
I ■ JR - IS ♦ 1 
DO 100 la IS. J 

BR(I) a BR(I) - BK(I*K)*D 
100 COHTXHUE 
110 JO a JR 
GOTO 90 
C 

C WARNING FORMATS 

2000 FORMATC/M ! WARNING ! ! 1 - IS SOLVER, LOSS OF AT LEAST 7 DIGITS' 

♦ /18X, 'Of REDOC Hi G DIAGONAL OF EQUATXQH; ' ,4X,I5) 

2001 FORMAT(/M! WARMING !! 2 - XH SOLVER, SIGN OF DIAGONAL CHANGED* 

♦ /18X, 'WHEN REDUCING EQUATION; ' ,181,15) 

2002 FORMAT (/' ! ! WARNING 1 J 3 - IN SOLVER. REDUCED DIAGONAL IS EERO' 

♦ /18X, 'FOR EQUATION; *,2SX, IS) 

2003 FORMATC/M! WARNING M 4 - IN SOLVER, RANK FAILURE FOR A ZERO' 

♦ /18X, 'UNREDUCED DIAGONAL IN EQUATION; ' ,7X, IS) 

2004 FORMATC/M! WARNING ! ! S - IN SOLVER, ZERO LOAD VECTOR') 

C 

END 

C=End Fortran 
C=Modul« DATEST 
CaBlock Fortran 

C SUBROUTINE DATEST(A, JH.DAVAL) 

recursive SUBROUTINE DATEST(A, JH.DAVAL) 


c+ 



*C 

c 



C 

c 

TEST FOR RANK 


C 

c 



C 

c 

INPUTS; 


C 

c 

A(J) 

- COLUMN OF UNREDUCED ELEMENTS IN ARRAY 

C 

c 

JH 

- NUMBER. OF ELEMENTS IN COLUMN 

C 

c 



C 

c 

OUTPUTS; 


C 

c 

DAVAL 

- SUM OF ABSOLUTE VALUES 

C 

c 

c+ 



C 

+C 


C 

C ARGUMENTS 

REAL*8 A(l) , DAVAL 
INTEGER JH 

C LOCAL ARGUMENTS 

INTEGER J 
C 

DAVAL > O.ODO 
DO 10 J a i, JH « 

DAVAL a DAVAL ♦ DABS(A(J)) 
10 CONTINUE 
C 

RETURN 

END 
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c 

c 

Cs£od Fortran 

OKodulc 007 

C*81ock Fortran 
C FUNCTION D0T(A,B,H) 

recursive FUHCTIOH DOT(A.B.H) 


C* r — +C 

C PURPOSE: C 
C THIS FUHCTIOH SUBROUTXHE PERFORMS THE DOT PRODUCT OF TWO C 
C VECTORS. C 
C C 
C ARGUMENTS: C 
C A - FIRST VECTOR IHVOLVED IH DOT PRODUCT C 
C B SECOHD VECTOR IHVOLVED IH DOT PRODUCT C 
C H HUMBER OF ELEMENTS IH EACH OF THE TWO VECTORS C 


am*8 dot, r(a) 

INTEGER H 


DOT ■ 0.0 
DO 10 I » 1, H 

DOT ■ DOT ♦ A(I)*B(X) 
10 COHTIHUE 
C 

RETURN 

END 

CaEnd Fortran 


File: nophlag.f 


CaModalo HOPHLAG 
CaAuthor K. Alvin 
CaOata July 1990 
CsBlock Fortran 

C 

C Subrout in* HOPHLAG 

C 

C Purpoa* : 

C This subroutine solves for the structural displacement 

C and valocity vactora at tha hali-step lor the phase lag 

C correction loop, and gate nev maasuraaant 

C 

C 

C Arguments 

C dalbata - dalta * bdamp 

C 

subroutine HOPHLAG (zp) 

raal *8 zp(l) 
include ’ shared . inc ’ 
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c 


LOCAL VARIABLES 


iatagar i 

saal*« v(MAXDOF) .dalbota 
C LOGIC 

C ABB APPLIES FORCES TO RHS AMD PREPARE MASS MULTIPLIER 

da tO i>t,adof 

gs(i) ■ gs(i) ♦ f (i) 

v(i) » (t. * dalta*ada«p ) *q( i ) * dalta«qdot(i) 

10 eoatiaaa 

C SOLVE FOR RIGHT HARD SIDE, ga 

do 77 i » i,adoi 

gs(i) * d«lsq*ga(i) ♦ v(i)*masa( jdiag(i)) 

77 cm&xsatM 

If (bdamp .n« . 0.) t&aa 
dalbata ■ dalta«bdasp 

call PHVKAO(stif , jdiag.q.ndof , dalbata, gs , 1 .dO) 

•adif 

C SOLVE FOR DISPLACEMENT, q, OS IMG RHS AMD MATRIX S 

f 

call SOLVER(as,ga,jdiag,adof ,2) 

do 100 i*l,ndof 

»(i) » (ga(i)-q(i))/dalta 
100 coatiaaa 

call ZEROVECT(zp,aa«a) 

do 200 jj » l.hdval 
i » hdrow(jj) 
j * bdcol(jj) 

zp(i) « zp(i) ♦ bd(jj)*ga(j) 

200 coatiaaa 

do 2S0 jj « l.hmral 
i » baroa(jj) 
j » hacolCjj) 

zp(i) ■ zp(i) ♦ bv(jj)*v(i) 

2S0 coatiaaa 

ratara 

aad » 


File: . zerovect . f 


CsModala ZEROVECT 

CaParposa Iaitializa vactor of givaa loagtb to zaro 
Calathor X. Alvia 
CaQata May 1990 
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* r 


C 

C Subroutine ZEROVECT 

C 



C 

subroutine ZEROVECT(v.n) 

real*8 v(l) 
integer n 

do 100 isl ( n 
r(i) a o.dO 
100 continue 

return 

end 


File: lu.f 


SUB ROOT lire LUFACT(A,N, PIVOT .DET.IER, UMAX) 
0****«*e**«****«****«»«**M****«****»****»**«*«»«*******«»*«*»*< 

C 

C SOBROOTIHE FACTOR USES GAOS SI AH ELIMIHATION WITH 

C PARTIAL PIVOTIHG AHD IMPLICIT SCALING TO DETERMINE 

C THE L*U DEC0MP0SI0H OF A SQUARE HATRIX "A" OF 

C ORDER N. THE ALGORITHM ALSO FINDS THE DETERMINENT 

C OF "A". UPON COMPLETION, THE ELEMENTS OF THE UPPER 

C TRIANGULAR MATRIX "U" ARE CONTAINED IN THEIR RESPECTIVE 

C LOCATIONS IN MATRIX "A". THE ELEMENTS OF HATRIX "L" 

C ARE CONTAINED IN THE LOVER TRIANGULAR PORTION OF "A", 

C BUT ARE SCRAMBLED WITH RESPECT TO "U" BECAUSE OF ROW 

C INTERCHANGE OPERATIONS NOT PERFORMED ON THE ELEMENTS 

C OF "L". THE VECTOR PIVOT (SEE BELOV) MUST BE USED TO 

C UNSCRAMBLE "L" IF IT IS TO BE USED FOR OTHER OPERATIONS. 

C 

C VARIABLES: 

C 
C 
C 
C 
C 
C 

c 
c 
c 
c 
c 

C*««*****«**«**e*e 

INTEGER PIVOT(l) ,IER,N,I, J,K,IO,NMAX 
REAL*8 A(NHAX,1) ,S(1000) ,C(1000) .DET.TEMP 
DET«1.0D0 

CM***W«SM*M**M»**M»*»**m»**«»***M**«**»*****MWM***M* 

C 

C FIND THE ROV NORMALIZING COEFFICIENTS S(I) FOR IMPLICIT SCALING. 

C EXIT ROUTINE IF ANY S(I)=0.0 


A=FULL SQUARE MATRIX (DOUBLE PRECISION) 

REORDER OF MATRIX A (INTEGER) 

PIVOT= VECTOR CONTAINING A RECORD OF 
ROV INTERCHANGES. THE INTEGER 
VALUE PIVOT(X) IS THE ROV VHICH 
VAS INTERCHANGED VITH ROV X AT 
FORVARD ELIMINATION STEP K. (INTEGER) 
DET»DETERHINENT OF MATRIX A (DOUBLE PRECISION) 
IER=ERROR FUG. IF IER=1, THE MATRIX A VAS FOUND 

TO BE SINGULAR, AND THE ROUTINE VAS EXITED. . IF 
IER=0, THE DECOMPOSITION VAS SUCCESSFUL. 
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nano 


C*®*®®®®*****^,* ®®«®®*®*®«*»tt®®*«®e®«*«®c»**«®a®s®®***®*®*®<!ii!<sa»*:S! 

DO 100 I»i,H 
S(I)»O.ODO 
DO 110 J*1,N 

IF (ABS(A(I,J)) .GT.S(I)) S(I)«ABS<A(I,J)) 

110 CONTINUE 

IF (S(I).Eq.O) THEH 
m«i 
DITsO.ODO 

asms 

END IF 
100 CONTINUE 

C®®**®®**®*®*®®*®*®*®®®®®®®*®®**®®**®®®*®**®*®**®®**®**®®*****®®* 


START FORWARD ELIMINATION STEP X 


i*®***«****®®®«****®® 

DO 130 X«l,Jhl 

C*®®®**®®®®*®®®®®*®®®*®®®**®*®®***®®®*®®*®*®®****®®*®®®®*®®®®®®®® 

c 

C DETERMINE PIVOT ELEMENT A(IO,K) BY FINDING THE ROW 10 

C BETWEEN K USD N CONTAINING THE MAXIMUM NORMALIZED . 

C VALUE IN COLUMN X. SET PIVOT(X)®IO 

C 

C®*®®®®®*®******®®®®®**®®®®®®***®®***®*®**®*®®**®*®®®®®®®*®®®®®®® 

C(K)«O.ODO 
DO 130 I*K,H 

TEMP«ABS(A(I,K)/S(I)) 

IF (TEMP. GT. C(lO) THEH 
C(K)»TEMP 
IO»I 
END IF 
130 CONTINUE 
PIVOT(K)®IO 

C®®*®®®®®®®®®®®®*®®®®®®®®*®®®®®®®®®*®**®®®*®®*®®®®®®®®®®®®*®®*®*® 

C 

C EXIT ROUTINE IF ALL VALUES IN COLUMN K AT OR BELOW 

C THE MAIN DIAGONAL ARE EqUAL TO 0.0 

C 

C*®®**®***®**®®®®®*®®*®®*®*®*****®®*®®®®®*®**®®®**®®*®*®®®®****®* 

IF (C(X).EQ.O.O) THEH 
IER»1 
DETaO.ODO 
RETURN 
END IF 

C**®**®®®®®**®®*®*®®®®®®*®®®®*****®®®®®®*®®*®*®*®®®®®®®®*®®®®®®*® 

c 

C INTERCHANGE ROWS 10 AND X FOR COLUMNS X TO N. SKIP IF IO=K. 

C SET DET«-DET IF ROWS ARE INTERCHANGED. 

C 

C**®**®*®***®®®*®*®®*®®®*®®®**®®®®®***®®*®®*®®®®***®®®*®*®®®®*®*® 

IF (lO.Eq.K) GOTO 150 
DET®-1.0D0*DET 
DO 140 J»K,N 
TEMP®A(K, J) 

A(K, J)®A(IO, J) 

A(IO,J)=TEHP 
140 CONTINUE 
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C*******************«*«*««******************v***«***!«e4*i3** «««:«** 

c 

C ELIMINATE counts X BELOV MAIN DIAGONAL BY MULTIPLYING 
C ROW K PROM COLUMH K TO H BY A(I ,K)/A(K ,K) AND SUBTRACTING 

C FROM ROW I. STORE THE MULTIPLIER FOR ROW I IN THE ELIMINATED 

C COLUMI X. MULTIPLY THE RUNNING PRODUCT DET BY DIAGONAL ELEMENT A(K,K). 

C 

C*************«*******«*a**»*****************************'**** , * , «** 


ISO DO 160 lag+l.N 

A(I,K)«A(I,K)/A(K,K) 

DO 170 J»K*1,N 

A(I,J)oA(I,J)-A(I,K)*A(K,J) 

170 CONTINUE 

160 CONTINUE 

DET»DET*A(K,K) 

120 CONTINUE 

C***********************.**************************************** 

c 

C CHECK LAST ROB/COLUMN FOR SINGULARITY. IF THERE IS NO ERROR, 

C COMPLETE CALCULATION OF THE DETERMINES! , SET THE ERROR FLAG 

C TO INDICATE NORMAL COMPLETION, AND EXIT. 

C 

C*****«***»*«******«*«*********«******************************«** 

IF (A(N.N).Eq.O.O) THEN 
IER»1 
DET=O.ODO 
RETURN 
END IF 

DET»DET*A(H,N) 

IER=0 

RETURN 

END 

SUBROUTINE LUSOLVE(A,N,B, PIVOT, NMAX) 
C**************************************************************** 

c 

C SUBROUTINE SOLVE 

C 

C**************************************************************** 
INTEGER PIVOT(l) ,N,I, J,K,NHAX 
REAL*8 A(NMAX,1), B(l), TEMP 
DO 100 Kal.N-t 

IF (PIVOT(K).EQ.K) GOTO 110 
TEHPaB(K) 

B(K)»B(J) 

B(J)aTEMP 

110 DO 120 I=K*1,N 

BCI)-B(I)-A(I,K)*B(K) . 

120 CONTINUE 

100 CONTINUE 

B(N)aB(N)/A(H,N) 

DO 130 IoN-1,1,-1 
DO 140 J»I*1,N 

B(I)»B(I)-A(I,J)*B(J) 

140 CONTINUE 

B(I)=B(I)/A(I,I) 

130 CONTINUE 
RETURN 
END 
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CsEHD Foams 
csBECI FACTA 

esPOaPQSg - Factors the A aatriz u L C « A, vita partial pivoting 
e*ASTSBE 8 X BELT IX, Sopt. 24 0 1987 
e 

e Input 

c aaat--— [ja X n] aatriz to bo factored, dastroyad on output 

e up — .-x—probla* siso 

e 

e Output 

e gnat—- contains tha LU dacoaposition 

c 

subroutino FACTA (anat ,np,arov,lp) 

c 

real*8 usat(*),eta 
tat agar lp(*) 

c 

do SO i«l,np 
SO lp(i)«i 
c 

c Find largest pivot . 
c 

do 100 k»l,np-l 
aaazkaO. 

BBSaxaO 

do 200 mek.np 

if (abs(aBat(lp(a)+(k-l)*nrou)) .gt. aaaxk) then 
aaaxks&bs ( anat ( lp(a) +(k-l ) *nrow) ) 


maz*s 

aadif 

200 continue 
l-lp(k) 

lp(k)alp(maaz) 

lp(anax)*l 

do 400 i=k+l ,np 

ata»anat ( lp( i ) ♦ (k- 1 ) *nrov ) /amat ( lp (k) ♦ (k- 1 ) »nr ow ) 
aaat ( lp( i ) •♦■ (k- 1 ) enrou) ■ et a 
do S00 jak+l.np 

anat(lp(i)*(j-l)*arou)saBat(lp(i)+(j-l)*nrou)- 
eta*amat(lp(k)+( j~l)*nrov) 

500 continue 

400 continue 

100 continue 


return 

and 


C EHD FORTRAB 
c=DECK LUSOLV 

c=PTJRPOSE - Solve L U x • b, 
c=AUTH0R 8 X' BELYIS, Sept. 24, 1987 
c * 
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c First solvss L j » b, than 0 x « y 
c 

c Input 

e mt*— [i X a] matrix factored by FACTA 

c np - — — problem six* 

c lp —pointer vector based on pivoting 

e rhs— RHS of aquation 

c 

c Output 

c aaat contains tha LU decomposition . 

c x— — — " the solution vector 

c 

subroutine LUSOLV(aaat ,np ,nrow , lp , rhs ,x) 
c 

raal*8 anat(e) ,x(*) ,rhs(») . 
integer lp(«) 
c 

do 54 i»l ,np 
x(i)<>rhs(i) 

50 continue 

c A - 

c Solve Lover System— — — — — — — — — 

c**a* Outer loop 

do 100 k*l,np 
rhs 00 =x (lp 00) 
if (rhs(k) . eq. 0.) go to 100 
c 

c**m Inner loop 
c 

do 200 i=k+l,np 

x(lp(i))=x(lp(i))-amat(lp(i)a(k-l)*nrov)*rhs(k) 
200 continue 
100 continue 
c 

c Solve Upper System— —————— — — — 

c 

do 300 k®np,l,-l 

x(k)=rhs(k)/amat(lp(k)+(k-l)anros) 
do 400 ial.k-1 

rhs(i)arhs(i)-x(k)*a«Bat(lp(i)+(k-l)*aroB) 

400 continue 

300 continue 
return 
end . 


File: prepcon.f 


C=Hodule PREPC0H 

C=PurpoBe Preprocess control analysis module for ACSIS 
C= Author K. Alvin 
C=Date June 1990 
CaBlock Fortran 

C 

C Subroutine PREPCOH 
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c 

C Purpose: ' 

6 This subroutine prepares ec, the control prediction integration 

C Matrix and eo, the observer construct matrix S (M+delta«‘D+delsq*K) 

C 

C 

C 

subroutine PREPCQS 

include ’ shared, inc' 

C LOCAL VARIABLES 

real* 8 me ,kc 
integer ier 

e Losie 

call READCQH 

C ! Form Control Prediction Integration Coefficient Matrix 
C 

C contype * -1 : Full State Feedback 

C contype » 0 Lnenberger Observer 

C contype a *1 Kalman Filter 

do 10 i ■ l.nact * naen 
do 20 j ■ l.nact ♦ nsen 
ec(i, j) ■ O.dO 
20 continue 

10 continue 

if (contype) 100,200,400 

J 00 ncsi * nact 

do 110 i ■ l.nact 
do 120 jj ■ l.bval 
j • bcol(jj) 
k » brov(jj) 

ec(i,j) » ec(i,j) ♦ delta*f2(i,k)*b(j j)/mass(jdiag(k)) 

120 continue 

110 continue 

goto 600 

200 ncsi * nact ♦ nsen 

do 210 i > l.nact 
do 220 jj » l.bval 
j * bcol(jj) 
k » brov(jj) 

•c(i, j) » ec(i.j) ♦ delta*f2(i ,k)*b( j j)/mass( jdiag(k)) 

220 continue 

do 240 j ■ aact+l.ncsi 
do 250 k > l.ndof 

ec(i.j) » ee(i.j) ♦ delta*f2(i,k)*12(k,j-nact) 
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250 

240 

210 


270 


290 

260 


400 


480 


S10 

800 

470 


830 


880 

520 

600 


1100 

C 


continue 

continue 

continue 

do 260 ii ■ l.hwal 
i ■ nact ♦ hvxow(ii) 
do 270 jj ■ l.bval 
j • bcolCjj) 
k » brou(jj) 

il (hvcol(ii) .ne. k) goto 270 

ec(i,j) ■ ec(i, j) ♦ delta*hv(ii)*b(jj)/naaa(jdiag(k)) 
continue 

do 290 j » nact*l,ncai 
k ■ kvcol(ii) 

ec(i,j) « ec(i, j) ♦ delta*h»<ii)*12(k, j-nact) 
continue 
continue 

goto 600 

ncsi s nact ♦ naen 

do 470 i 3 l.aact 
do 480 jj a l.bval 
j « bcol(jj) 
k » broe(jj) 

ec(i, j) a ec(i,j) ♦ deltaaf 2(i ,k)*b(j j)/maaa( jdiag(k)) 
continue 

do 800 j a nact+l,ncai 
do 510 k a i.ndof 

ec(i,j) a ec(i, j) ♦ delta*f2(i,k)»12(k, j-aact)/ 
mass( jdiag(k) ) 

continue 

continue 

continue 

do 520 ii a l.bwal 
i a nact + hvrou(ii) 
do 630 jj a i.bval 
j » bcol(jj) 
k » brow(jj) 

il (ivcol(ii) .ne. k) goto 530 

ec(i,j) » «c(i.j) ♦ dolta»hv(ii)*b( j j)/nass( jdiag(k)) 
continue 

do 850 j a nact>l,ncai 
k a bvcol(ii) 

ec(i,j) a ec(i, j) ♦ delta*hv(ii)*12(k, j-oact)/ 
naaa( jdiag(k)) 

continue 

continue 

continue 

do 1100 i * l.ncai 

ec(i,i) a «c(i,i) + l.dO 
continue 

FACTORIZE ec 

call FACTA(ec,ncai,KAXCSX .pivot) 


85 



if Cits . «q. 1) than 

prist « , ’ PREPCON : Singular Matrix for Control Integration* 
•ffidii 

C Fora 9b# error Integration Coefficient Matrix, eo 

k « 1. ♦ delta«adasp 
ke » delta«bdamp ♦ delsq 
do 1200 i»l.mlen 

ee(i) ® ac*aass(i) * kc*stif(i) 

1200 continne 

C FACTORIZE eo 

eall S0LYER(eo,go, jdiag.ndof ,1) 

C Initialize Observer States 

sail S8OTZCT (qe , sdef) 
eall ZEROVECT(qedot.ndof) 
eall ZEROVECT(pe.adof) 

return 

end 


File: control. f 


CsHodule CONTROL 
CsAnthor K. Alvin 
C=Date Hay 1990 
CaBlock Fortran 
C 
C 
C 
C 
C 
C 
C 
C 
C 
C 
C 
C 
C 
C 
C 


Subroutine CONTROL 
Purpose: 

This subroutine carries out the numeric integration of one time 
step of the control system. 


Arguments . 

q*p 

qedotp - 

PP 

z 


estimated displacement vector at half time step 
estimated velocity vector at half time step 
generalized momentum (f-D*qedotp~K*qep) 
measured sensor output 


subroutine CONTRQL(z) 


include * shared. inc* 
real*8 z(l) 

C LOCAL VARIABLES 


real*8 qep(HAXDOF) , qedotp (HAIDOF) .pp(HAXDOF) . v(HAIDOF) 
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C LOGIC 

C Fora BBS of Control Prediction Equation Sat 

C 

C contype ■ -1 Full State Feedback 

C contype ■ 0 : Lnenberger Observer uith Li = 0 

C contype ■ ♦! : Kalman Filter v/generalized momentum variable 

4 

call ZZBQVECKgc.ncsi) 

if (contype) 100,200,300 

100 continue 

do 110 i ■ l,ndof 

q*(i) ■ q(i) 

qedot(i) » qdot(i) 

110 continue 

200 , continue 

do 210 i ■ l.ndof 

qep(i) ■ qe(i) ♦ delta*qedot(i) 
qedotp(i) * qedot(i) 

pp(i) « f(i) - mass( jdiag(i))*adamp*qedotp(i) 
v(i) * qep(i) bdaop*qedotp(i) 

210 continue 

call PNVHA0(stif , jdiag, v.ndof ,-l .d0 ,pp,l .d0) 

do 220 isl.nact 
do 230 j » l.ndof 

gc(i) ■ gc(i) - fl(i,j)*qep(j) - f2(i, j)*(qedot( j) ♦ 
delta*pp ( j ) /mas e ( j diag ( j ) ) ) 

230 continue 

220 continue 

if (ncai . eq. nact) goto 600 

do 240 ianact+l.ncai 
k a i - nact 
gc(i) » z(k) 

240 continue 

do 245 ii » l,hdval 
i a bdrou(ii) ♦ nact 
j a hdcol(ii) 

gc(i) ■ gc(i) - bd(ii)*qep( j) 

24S continue 

do 250 ii a l'hwal 
i a kvrou(ii) ♦ nact 
j a hvcol(ii) 

gc(i) ■ gc(i) - hv(ii)*(qedot( j)>delta*pp( j)/masa( jdiag( j))) 
250 continue 

goto 600 

300 continue 

do 310 i a l.ndof 
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Vp(i) • q«(i) 

pp(i) » f (i) - aaBB(jdiagd))*adanp*qapd)/dalta 
▼(i) ■ (1 ♦ bdamp/dalta)*qapd) 

310 coatiaaa 

call P*7HJU>(«*if,jdiag,v,adof,-l.dO,pp,i.dO) 

do 320 i»l,aact 
do 330 j » l,adof 

gc(i) « gc(i) - 21(i,j)*q«p(j) - 22(i, j)*(pa(j) ♦ 

. dolta*pp(j))/Bass(jdiag(j}) 

330 coatiaaa 

320 eoatiaaa 

do 340 isaaet*l,acsi 
1 « i “ aaet 
ge(i) « *(k) 

340 coat lane 

do 346 ii ■ l.hdval 
i » bdxwdi) aact 
j « idcol(ii) 

gc(i) ■ ge(i) - bddi)*q«p(j) 

345 coatiauo 

do 350 ii ■ 1 ,hvval 
i « hvrop(ii) *■ aact 
j ■ hveol(ii) 

gc(i) ■ gc(i) - hv(ii)*(pa(j) ♦ dalta*pp(j))/ma.Ba(jdiag(j)) 
350 coatiaaa 

C FIHD r, C0HTR0L AHS STATE CORRECTION FORCES 

600 caU LOSOLV(ac,acai.HAXCSX. pivot, gc.r) 

do 610 j«l,aact 
a(j) • r< j) 

610 coatiaaa 

do 620 j*nact+l ,ncsi 
gaasa( j-aact) ■ r(j) 

620 coatiau* 

C FIHD CONTROL CONTRIBUTION TO RHS VECTOR FOR 
C OBSERVER AHD STROCTRE 

do 710 ial.adof 
gt(i) « O.dO 

god) ■ O.dO 

gk(i) * O.dO 

710 coatiaaa 

do 720 jjal.bval 
i » brow(jj) 
j » bcol(jj) 

gt(i) a gs(i) ♦ b(jj)*a<j) 

go(i) a go(i) ♦ b(jj)*a(j) 

gk(i) a gk(i) ♦ b(j j)*a(j) 

720 coatiaaa 

ii (coatjpa .aq. 0) than 
do 725 i > i.adof 

do 730 jai,asaa 

god) » god) ♦ maBB(jdiagd))*I2d,j)*gajnma.(j) 

730 coatiaaa 
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725 continue 

elseif C coat 7 p • .eq. 1) than 
do 73S i » l.ndof 
do 740 j»l,naen 

go(i) ■ go(i) + (12<i,j)anaaa(jdiagCi))aii(i,j) /delta) 
*gaaaa(j) 

gk<i) ■ gk(i) ♦ 12(i, j)*gaama(j) 

740 continue 

735 continue 

end if 

return 

end 


File: secorder.f 


CaModnle SECONDER 
C* Author K . Alvin 
ODate Nay 1990 
CsBlock Fortran 


C 

C 

C 

C 

C 

C 

C 

C 

C 

C 

C 

C 

C 

C 

C 

C 

C 

C 

C 

C 

C 

C 

C 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 


Subroutine SECORDER 
Purpoae: 

Solvea the second-order dynamical equation: 

Nr" ♦ Ox' ♦ b • 1 ♦ j 

at time (a+1) given f(n*l/2), g(n+l/2) and x,x' at n by 
the midpoint implicit integration rule. 

Step size is 2*delta. 

0 is of the form (alpha«H ♦ beta*K), f is an applied force, 
and g is assumed to be other applied force from a feedback 
control loop. The matrix E is the factored fora of the 
integration coefficient matrix: E a (X ♦ delta*D ♦ delta'2*K) . 


Arguments : 

m 

k 

alpha 

beta 

f 

S 

e 

x 

xd 

delta 

delaq 

jdiag 

ndof 


matrix N 

matrix K 

scalar alpha 

scalar beta 

Force vector fCn+1/2) 

Feedback force vector g(n*l/2) 
martix S 

Variable vector x(n) 

Variable vector x' (n) 

Half of integration time step 
delta‘2 

Diagonal location pointer for X,K,E matrices 
Humber of equations and length of q 
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C t mass multiplier for RHS preparation 

C delbeta - delta * beta 

C subroutine SECORDER(a, k, alpha, beta, i,g,e,x,xd, 

C delta, delsq,jdiag,adof .KAX9QF) 

recursive subroutine SECORDER(a,k, alpha, beta.f ,g,e,x,xd, 
delta, delsq, jdiag.ndof .MAXDOF) 

C ARGUHESTS 

real *8 a(l),k(l), alpha, beta, f (1} ,g(l),e(l) 
real«8 x(l) ,xd(l) , delta, delsq 
integer jdiag(l) ,ndof .HAIDOF 

C LOCAL VARIABLES 

integer i 

real *8 ▼ (3000) ^delbeta 

C LOGIC 

C ADD APPLIED FORCES TO KHS ADD PREPARE MASS MULTIPLIER 

do 10 i»»l,ndoi 

g(i) ■ g(i) + f(i) 

v(i) ■ (1. ♦ delta«alpha)«x(i) ♦ delta*xd(i) 

10 continue 

C SOLVE FOR RIGHT HAHD SIDE, g 

do 77 i»l ,ndof 

g(i) a delsq*g(i) ♦ v(i)*m(jdiag(i)) 

77 continue 

if (beta .ne. 0.) then 
delbeta » delta*beta 

C Activate EBE computations lor internal force by using STIFFRC 
C subroutine. Otherwise use PKVHAD (profile matrix/vector mult-add 

C call PHVHAD(k, jdiag,x,ndof .delbeta, g, 1 .dO) 

call STIFFRC(x, delbeta, g) 
endif 

C SOLVE FOR DISPLACEHEHT, q, US IMG RHS AMD MATRIX E 
call SQLVER(e,g, jdiag.ndof ,2) 
do 100 i>l,ndof 

xd(i) ■ 2.*(g(i) - x(i))/delta - xd(i) 
x(i) ■ 2.*g(i) - x(i) 

100 continue 

return 0 

end 


File: measure. f 
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CaHodule MEASURE 
CaAuthor K. Alvin 
C=Data Kay 1990 
Csfllock Fortran 

C 

C 

C Subroutine MEASURE 

C 

C Purpose : 

C This subroutine storoa now measured aanaor data by using the 

C previous displacement and velocity vectors at the sensor 

C locations 

C 

C 

C Arguments 

C zp - measured sensor data array 

C 

subroutine KEASURE(zp) 

include ’ shared . inc ' 
real*8 zp(l) 

call ZERQVECT(zp,nsen) 

do 100 jj » l.hdval 
i » hdrow(jj) 
j ■ hdcolCj j) 

zp(i) ■ zp(i) ♦ hd(jj)*q(j) 

100 continue 

do 200 jj ■ l.hwal 
i » hvrow(jj) 
j * hvcol(jj) 

zp(i) » zpCi) ♦ hv(jj)*qdot(j) 

200 continue 

return 

end 


File: eigens.f 


C=Module EXCESS 

C=Purpose Find Eigenmodes given Mass, Stiffness Matrices 
C=Author K. Alvin 
CsDate March 1990 
C=Block Fortran 

C 

subroutine EXCESS 
C 

C 

C 

C COKHOS ABD GLOBALS 
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« 4 


iaeloda ’sharad.iac' 

C LOCAL VARIABLES 

para*atar(NVI*100 ,HHNVH*NAXDOF*HVM ,HXCHVsNVH*(NVM*l)/2) 
iatagar ialCNAXDOF^awae.i.j.k.kk^oat 
?«al*8 vl(NNNVN) ,vx(HNNVH) ,akk(NXCNV) t &ffl»(ilXCNY) 
saal*8 xx(NVH,NVl) ,eigv(NVH) ,aigold(NVH) 
real*8 tolaig,toljae, omega, fhz 

toleiga.OQOl 
owt * 13 

toljae ■ toleig 
avoe * mia0(2*aaig,100) 
awac « miaO(avec,adof ) 

C SET Iff ISL VECTOR 

iel(l) a i 
do 50 j«2,sdo£ 

i»Kj) a j - jdiag(j) ♦ jdiag(j-l) *1 
50 continue 

C CALL EIGENSOLVER 

call SSPACECstii.maaa.vl.irr.akk.aaffls.zr.eigv.aigold.isl, 
j diag , aaig ,av«c , ndof . tolaig , tol j ac , out ) 

C WRITE OUTPUT 

writaCoat.*) *EIGEH ANALYSIS RESULTS:* 

writaCout,*) • RADIAL CYCLIC* 

writ a (out, *) * MODE EIGENVALUE FREQUENCY FREQUENCY* 

writaCout ,*) * — — - — — --* 

writa(out,*) 
do 100 ial.aaig 

omaga a dsqrt(eigv(i) ) 

ihz a omaga/(2*3. 141692654) 

writaCout,’ (i5,3(3x,gl2.5))’) i,eigv(i) .omega, Ihz 
100 coatiaaa 

writ«(out,*) * EIGENVECTORS : ’ 
do 200 j*l,naig,5 
writaCout,*) 
do 300 i*l,ndoY 

k. « ndof*(j-l) ♦ i 

writaCout, *(i5,5(lx,gl2. 5))') i,(irr(kk) ,kk=k,k+4*ndof ,ndof) 
300 coatiaaa 

200 coatiaaa 

writaCoat,*) 'NASS MATRIX DIAGONAL:* ' 
do 400 i*l,nnp 
do 450 j*l,6 
if (id(j ,i) .ne.O) than 

writaCoat,*) i, j ,id(j ,i) ,mae«(jdiag(id(j ,i))) 
endif 

450 coatiaaa 
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400 


continue 


return 

end 

C=End Fortran 


File: singaig.f 


SUBROUTIHE SSPACE (AK,AM,VL,VR,AKK,AMM,XX,EIGV,EIGCLD ,ISL, 
1 IDIAG .HEIG.NVEC , HD0F , TOLEIG , TOLJ AC , HV) 


c 

e 

e 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

e 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

e 

c 

c 

c 

c 

c 

c 


input : 


AK 

AM 

TSL 


IDIAG 

MEIG ; 

HVEC 

NBOF 

TOLEIG 

TOLJAC 

HV 


stiffness matrix (profile values) (NDOF) 

consistent mass matrix (profile values) (HDOF) 
stores in position “i" the -row # of tip of 
column "i" (NBOF) 

position of diagonal terms in profile (HDOF) 

S of required eigenvalues 

9 of subspace vectors 

9 of degrees of freedom 

tolerance for eigenvalues convergence 

tolerance for Jacobi convergence 

logical unit number for output 


output : 

VL(HDOF.HVEC) 

- VL(.,1..MRH0D) : 

- VL(. .HRMQD. .HVEC) : 

VR(HDOF,HVEC) : 

- VR(.,1.. HRMQD) : 

- VR(.,HRM0D'. .HVEC) : 

AXX( HVEC* (HVEC ♦ l)/2) 
AMM(HVEC*(NVEC ♦ l)/2) 
XX(HVEC*HVEC) 

EIGV(HVEC) 

- EIGV(1..HRM0D) : 

- EIGV(HRM0D . .HVEC) : 


working array 
AM times rigid modes 
subspace at the previous step 

eigen-vectors 
rigid modes 

subspacs at this step (eigenvectors) 

stiffness matrix in the subspace 
consistent mass matrix in the subspace 
subspace eigenvectors 

current eigenvalues 
0 eigen-values 
following eigenvalues > 0 


EIGOLD(HVEC) : same as EZGV 


IMPLICIT REAL *8 (A-H.O-Z) 

DIHEHSIOH AK(1) ,AM(1) .VL(HDOF.HVEC) ,VR(NDOF,HVEC) ,AKK(1) , 
1 AMX(i) , II (HVEC* HVEC) ,EIGV(1) ,EIG0LD(1) ,ISL(1) , IDIAG(l) 

C 

VRITE (HV.1003) HEIG .HVEC, HDOF, TOLEIG 
C 

CALL IHVECT ( AK, AM, VL.VR, IDIAG, HDOF, HVEC) 

CALL FACT (AK, IDIAG, ISL.HDOF.HV) 

CALL HULL (AK, AM, VL.VR, IDIAG, ISL.HDOF.HRMOD) 
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c 

c HHHQD : S of rigid sodas 

C 

KRXTE (STB, 1004) HRMOD 
ffSUl*«VEC-IRIOD 

C - 

c m ORTHO (VR.YL.HDOF.HRSOD.HVEC) 

C 

SXtsO 

SSSAXsjg 

®VECl»lVEC-i 
DO S I»1,SVEC 
6 RIGOLS(I)»0.0 
C 

500 KXT*Hrm 

WRITE (8W. 1000) HIT 

C 

CALL SOLVES (AK ,V 1U ^BUBBKl) .WRCUMUmD*! ) .J3X1G . ISL.HDOF^HSOB) 
CALL ORTHO (VL.TR, HDOF, HRMOD, HVEC) 

C 

XJ*0 

DO 10 J*HRK0D«-1,HVEC 
C 

C CALCULATE THE UPPER PART OP AXX (SYMMETRIC) 

C 

DO 10 IaNRMODal.J 
TR*0.0 

DO 11 K-l.HDOF 

11 TR»TR0VL(K,I)*VR(K.J) 

IJ=IJH 
AKK(IJ)aTR 
10 CONTIHUE 
C 

CALL MULT (AM.VRd.HRMOD+D.VLCl.HRMOD+D.ISL.IDIAO.NSUB.NDpF) 

C 

IJaO 

DO 20 JaHRMOD+1 ,NVEC 
C 

C CALCULATE THE UPPER PART OF AMM (SYMMETRIC) 

C 

DO 20 IaHRMOD+1 . J 
TRaO.O 

DO 21 K»1,HD0F 

21 TR“TR+VL(K,I)*VR(X, J) 

, IJaIJ+1 
AMM(IJ)*TR 
20 CQHTIHUE 
C 

CALL JACOBI (AKX,AHH,XX t EIGY(HRMOD+l),NSMAX,TQLJAC,SSUB,NV) 

C 

C ORDER EIGENVALUES k EIGEHVECTORS 
C 

30 IS-0 

DO 40 I»HRH0D*1,HYEC1 

IF (EIGY(I*1) .GE.EIGV(I)) GO TO 40 
IS«1 

TR«EIGV(I*1) 

EIGV(I*1)=EIGY(I) 
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EIGV(I)=TR 
DO 41 J-i.HSUB 

TBaXX(J*(I"HRMOD)*HSUB) 

XX(J*(I-HRMOD)*HSUB)aXI(J>(I-HMOD-l)*HSTO) 

41 XX(J*(I-HRM0D-1)*HSOT)»TS 

40 CONTINUE 

IF (IS.EQ.l) GO TO 30 
C 

C SUBSPACE CONVERGENCE TEST 
C 

IC0HV»0 

DO 60 I»HRHOD'H,HVEC 

TRaDABS((EIGOLD(I)-EIGV(I) )/EIGV(I) ) 

El GOLD ( I ) =EIGV (I) 

EIGV(I)*TR 

IF (TR.GT.TOLEIG.AHD.I.LE.NEIG) IC0HV»1 
SO COHTIHUE 

WRITE (HW.1001) (EIGV(I) , I»1,NVEC) 

if (laarv.Eq.o) go to ioo 

c 

IF (HIT . LE . NITNAZ) GO TO 70 
WRITE (NV.1002) 

GO TO 100 
C 

C UPDATE EIGEN VECTORS 
C 

70 DO 80 Ial.NDQF 

DO 80 J-l.HSUB 
TR*0 . 0 • 

DO 81 K=1,HSUB 

81 TR=TR+VR(I,K*NRH0D)*XX(K*(J-1)*HSUB) 

VL ( I , J ♦NRMOD ) =TR 
80 COHTIHUE 

DO 90 I*1,HD0F 

DO 90 J*NRM0D+1 ,HVEC 
90 VR(I,J)»VL(I,J) 

GO TO 600 
C 

C CALCULATE FIHAL EIGENVECTORS 
C 

100 DO 110 I*1.HD0F 

DO 110 J»1,HSUB 
TR*0.0 

DO 111 K*1,NSUB 

111 TR*TR* VL(I , K+NRMOD ) * XX ( K* ( J- 1 ) »NSUB ) 

VR(I,J+HRHOD)»TR 

110 COHTIHUE 

DO 112 I»1,NVEC 

112 EIGV(I)aEIGOLD(I) 

C 

RETURN 

C 

1000 FORHAT ( SX , 12HITERATI0H HO, 15) 

1001 FORMAT (6(2X,1PE10.3)) 

1002 FORMAT (5X.24HWE ACCEPT CURRENT VALUES) 

1003 FORHAT (//20X, 'SUBSPACE ITERATION ROUTINE'//' NB OF EIGENVALUES*' , 
1 IS/' SB OF VECTOR* ',16/* HB OF DOF* ',15/' TOLERANCE* ' , 1PE10 . 3/) 

1004 FORMAT (• HB OF RIGID MODES** .IS//) 

C 
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END 

cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc 

SUBROUTINE INVECT (AX, AM.VL, VR.IDIAG .NDOF.NVEC) 

IMPLICIT REAL*8 (A-H.O-Z) 

DIMENSION AI(i),,AM(l),VL<l),VR(NDOF,NVIC) e imG(i) 

C 

hd-ndof/hvec 

c 

DO 10 :>i,NDQF 

neisusd) 

ra(i,i)«AM(xx) 

n.(i)»AM(ii)/«(ii) 

DO 10 J-2.NVEC 
VR(I,J)a0.O 
10 ccnmsuE 
c 

LLaNDOF-ND 

C 

DO 20 4*2»JREe 
TR=0 . 0 
C 

DO 30 I»1,LL 

IF (VL(I).LT.TR) GO TO 30 

TRsVUX) 

IJ«I 

30 CONTINUE 

C 

DO 40 I»LL,NDOF 

IF (VUl) .LE.TB) GO TO 40 
TRaVL(I) 

IJ»I 

40 CONTINUE 

C 

VL(IJ)»0.0 

LL»LL-ND 

VR(IJ,J)»1.0 

C • v 

20 CONTINUE 
C 

RETURN 

END 

CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC cc 
SUBROUTINE MULT (AM,VR,VL,ISL,IDIAG,HVEC,NDOF) 

IMPLICIT REAL*8 (A-H.Q-Z) 

DIMENSION AM(l),VR(NDOF,HVEC),VL(NDOF,HVEC),ISL(l) f IDIAG(l) 

DO 500 IVal.HVEC 
DO 100 I»1,ND0F 
TR=>0.0 
IJalDLAG(I) 

IK»ABS(ISL(I>) 

KK»I 

DO 110 X>IX,I 

TR=TR+AM(IJ)*VL(KK,IV) 

IJaIJ-1 

KK=KK-1 

110 CONTINUE 

C 

IF (I.Eq.NDOF) GO TO 99 
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IX«H-1 

DO 120 K^IK.HDOF 

I? (I.LT.ABS(ISL(K))) CO TO 120 
IJ«IDIiG(K)-K4.I 
TR*TR+AM(IJ)*VL(K , IV) 

120 CONTINUE 

99 VR(I , IV)*TR 

100 CO N T INUE 

500 COHTIHUE 

HETURH 

END 

cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc 

SUBROUTINE JACOBI (AK,AH,XX,EIGV,BSHAX,TOL,N,NV) 

IHPLICIT REAL *8 (A-H.O-Z) 

DIMENSION AK(l),AM(l),XX(N,N).EIGV(l) 

C 

C INITIALIZE 
C 

ZEROoO.O 
TW0«2.0 
DO 10 I»1,H 

II*I*(I-1)/2*I 

IF (AK(II).LE. ZERO. OR. AK(II).LE. ZERO) GO TO 900 
EIGV(I)»AJC(II)/AM(II) 

DO 20 J°1 ,N 
XX(I,J)=ZERO 
20 CONTINUE 

XX(I.I)»1.0 
10 CONTINUE 
C 

C SET COUNTER 

C 

NSWEEP»0 

NR=N-1 

SOO NStfEEP*»NSVEEP+l 
C 

C CHECX IF ZEROOING IS REQUIRED 

C 

EPS=0 . 0 1**NSWEEP 
EPS=EPS*EPS 
DO 150 J»1,NR 
IIK»J*1 

DO 150 XallX.N 
JJ»J*(J-1)/2*J 
KK«K* (K- 1 ) /2*K 
JK=K*(K-l)/2+J 

EPSAK«( AK ( JK) * AK ( JK ) ) / ( AX ( J J ) * AK (KK) ) 
EPSAM«(AH(JK)*AM(JK))/(AH(JJ)*AM(KK)) 

IF (EPSAX.LT.EFS.AND.EPSAM.LT.EPS) GO TO 150 
C 

C CALCULATE ROTATION ELEMENTS 

C 

AKK»AK(KK)*AM( JK)-AM(KK)*AK(JK) 

AJJ»AK(JJ)*AX(JK)-AH(JJ)*AK(JK) 

AB=(AK(JJ)*AM(XX)-AX(XX)*AM(JJ))/TVO 

CHECK=AB*AB+AKX*AJJ 

IF (CHECK.LT. ZERO) GO TO 900 

SqCH=DSQRT( CHECK) 
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D1*AB+SQCH 

D2alB“SqCH 

DEH»D1 

IF (DABS(D2) .GT.DABS(Dl)) DEH=D2 
I? (S19.lS.ZZaO) GO TO 46 
CA«ZEBO 

CGs»AS(JI)/AX(KK) 

GO TO SO 

46 ClaAKX/DEH 

GGs-AJJ/DEH 


PMFtffil 6EHERALIZED BOTH I OH 


SO 


60 

70 


80 

90 


100 

110 


JPlaJ+t 

JM1»J-1 

KP1»K*1 

KH1»X~1 

IF (Jli.LT.l) GO TO 70 

00 60 

IJ s J*JMl/2+I 

IKaK*KM/2*I 

AKJaiK(IJ) 

AKX-AX(IK) 

AMJaAM(IJ) 

AXK»AM(XX) 

AX(IJ)=1XJ*CG*AXX 

AH(IJ)=AJU+CG*AHX 

AX(IK)=AXX*CA*AXJ 

AM(IK)»AHK«-CA*AJIJ 

C08TIBUE 

IF (KP1.GT.H) GO TO 90 
DO 80 I»XP1,H 
JI*»I*(I“l)/2+J 
KI»I*<I-1)/2*X 
AKJ«AK(JI) 

AMJaAM(JI) 

AKK 3 AK(KI) 

AKX 3 AH(XI) 

AK(JI) a AXJ+CG*AXX 

AH(JI) 3 AHJ*CG*AMK 

AX(XI) 3 AXX+CA*AXJ 

AM(XI) 3 AHK+CA*AHJ 

COHTXNUE 

IF (JP1.GT.XH1) GO TO 110 
DO 100 I-JP1.XH1 
JI«I*(I-1)/2*J 
IX*8>(X“1)/2*I 
AXJ»AX( JI) 

AMJ»AM(JI) 

AKK=AX(IX) 

AMKbAN(IK) 
AX(JI)«AXJ*CG*AXX 
AH( JI)»AXJ+CG*AKX 
AX(IK) 3 AKK+CA*AKJ 
AX(IK)*AHK+CA*AMJ 
CQBTIHUE 
AKX«AX(XX) 

IKK 3 AX (XX) 

AXJ a AX(JJ) 
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AXJaAM(JJ) 

AK(KX)»1KKVTV0*CA'*AX(JK)+CA*CA*AKJ 

AM(KX)«AXX+>TVO*CA*AH(JX)*CA*CA*AXJ 

AX ( J J ) a AX J ♦TVO*CG* AX ( JX ) •*-CG*CG« AKK 

AM(JJ)»AMJVnfO*CG*AM(JK)*CG*CG*AKK 

AX (JX) “ZERO 

AN(JK)»ZERO 

UPDATE EIGEHVECTOR FOR THIS ROTATION 


120 

150 


DO 120 I»1,N 
XXJ«XX(I,J) 
XXK a XX(l,K) 
XX(I t J)aXXJ*CG*XXX 
XX(I f K) a XXK*CA*XXJ 
CONTINUE 
CONTINUE 


UPDATE ggEHVALDES * CHECK CONVERGENCE 

NT«N*(N-H)/2 
ICOHVaO 
DO 160 I»1,H 

II*I*(I-l)/2+I 

IF (AK(II).LE. ZERO. OR. AH(II).LE. ZERO) GO TO 900 
TR=AX(II)/AM(II) 

DEH» (TR-EIGV ( I) ) /TR 
EIGV(I)*TR 

IF (DABS(DEN) .GT.TOL) IC0NV=1 
160 CONTINUE 

IF (ICOHV.Eq.l) GO TO 499 

CHECK OFF DIAGONAL TERNS 

EPS*TQL*TOL 
DO 170 J»1,NR 
IIKaJ+1 

DO 170 K»IIK.H 
JJaJ*(J-l)/2+J 
KX*K*(K-l)/2+K 
JK»K*(K-l)/2+J 

EPSAK a (AK(JK)*AK(JK))/(AK(JJ)*AK(KK) ) 
EPSAN»<AH(JK)*AH(JK))/(AM(JJ)«AX(KK)) 

IF (EPSAK.LT.EPS.AND.EPSAN.LT.EPS) GO TO 170 
GO TO 499 
170 CONTINUE 

SCALE EIGENVECTORS 


179 DO 180 I»1,N 

II»I*(I-l)/2+I 

AXK*DSQRT(AM(II)) 

DO 180 Jal.N 

XX(J,I)»XX(J,I)/AKK 

180 CONTINUE 

* 

RETURN 

* 

499 IF (NSVEEP.LE.NSNAX) GO TO 500 
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WRITE (NV.1000) 

GO TO 179 
C 

900 WRITE (IV, 1001) 

STOP 

C 

1000 FORMAT (6X,34BM0 CONVERGENCE AT BSMAX ITERATIONS) 

1001 FORMAT (SX.46HERR0R IB JACOBI : XATRII HOT POSITIVE DEFINITE) 



cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc 

SOiaOOTIBE SOLVES (AI.VL.VR.IDIAfi.ISL.BDOF.NVEC) 

C ? 

C PURPOSE : SOLVES THE SINGULAR PROBLEM : AX x VL ■ VR 
C TEE SINGULAR COLUMNS INTO AX ARE INDEXED BY 

C TEE NEGATIVE VALUES OF ISL, TEE CORRESPONDING 

C TERMS OF THE SOLUTION ARE PUT TO 0 

C 

c 

IMPLICIT REAL*8 (A-H.O-Z) 

DIMENSION AX(l) ,VL(NDOF,NVEC) ,VR(NDOF,NVEC) ,IDIAG(l) ,ISL(1) 

C 

DO SOO IV»1,HVEC 

DO 60 I»1,ND0F 
60 VL(I,IV)»VR(I,IV) 

BACKSUBSTITUTE 

DO 100 XC>2,HD0F 
TRsO.O 
ICl-IC-1 

IH1«IDIAG(IC)-IC 
IKalSL(IC) 

IF (IX.LE.O) THEN 
VL(IC,IV)aO. 

GOTO 100 
END IF 

IF (IK.GT.IC1) GO TO 100 
DO 120 X-IX.IC1 

TR*TR+AK( IN1+K)*VL(K , IV) 

CONTINUE 

VL(IC,IV)=>VL(IC,IV)-TR 
CONTINUE 

SOLVE DU»U 

DO ISO ICal.HDOF 

VL(IC.IV)3VL(IC,IV)/AK(IDIAG(IC)) 

ISO CONTINUE 

BAK SUBSTITUTE 

IICaBDOF 
DO 200 IC»2,HD0F 
TRaVL(IIC.IV) 

IC1»IIC-1 
IXalSL(IIC) 
nri»iDiAG(iic)-iic 
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IF ((IK.GT.IC1).0R.(IX.LE.05) GO TO 221 
DO 220 K*IK,IC1 

VL(K ,IV)=VL(K,rV)-lX(ISl+K)*TR * 

220 CONTINUE 

221 IICaIIC-1 
200 COBTXSUB 
BOO CONTINUE 

RETURN 

EHD 

cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc 

SUBROUTINE FACT (AX,IDIAG.ISL,NDOF,HV) 

c 

C PURPOSE : LDLT DECOMPOSITION OP THE POSITIVE SERI-DEFINITE 
C MATRIX AX, THE SINGULAR COLUMNS OF AX ARE INDEXED 

C BY A NEGATIVE VALUE OF IDIAG> 

C 

C 

UEELICTT B£AL«8 (A-H.0-Z3 
DIMENSION AX(1),2DIAG(1),ISL(1) 

C 

C DETERMINE MIN A MAX 
C 

trsdabs(akcD) 

AMIN-TR 

AMAX*TR 

DO 10 IL»1,ND0F 

TR»DABS(AK(IDIAG(IL))) 

IF (TR.LT.AMIN) AMIN“TR 
IF (TR.GT.AMAZ) AHAX»TR 
10 CONTINUE 

ZER0 3 (AMAZ+'AMIH)*1 .0D-10 

c 

C LOOP OVER COLUMN 
C 

DO 100 ICal.NDOF 
MIC=ISL(IC) 

IClaIC-1 

MIC1=HIC-*T 

IN2*IDIAG(IC)-IC 

IF ((HIC.LT. 1) .OR. (MIC.GT. IC) ) GO TO 901 
C 

C CALCULATE GS 

C 

IF (MIC1.GT.IC1) GO TO 150 
DO 120 ILoHICl'ICl 

IF (XDIAG(IL) .LT.O) THEN 
AX(IN2+IL) 3 0 
GOTO 120 
END IF 

MILalSL(IL) 

ILl=IL-i 

MIH»MAXO(MIL.MIC) 

IH1*IDIAG(IL)-IL 

TRaO.O 

IF (MIM.GT.IL1) GO TO 120 
DO 130 KaMIM.ILl 

TR=TR+AK(I51+K)»AXCIS2+K) 

130 CONTINUE 
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IN»XH2*IL 
AK(XM)*AK(IN)-TR 
120 CONTINUE 

C 

C CALCULATE LAO 

e 

ISO TR»0.0 

IF (IIC.6T.IC1) 60 TO 201 
DO 200 ILaMIC.ICl 

IF (IDIAG(IL).LT.O) GOTO 200 
AG*AX(XS2*IL) 

AL*AG/AX(XDIAG(IL)) 

AK(IH2*IL)«AL 

TR“TR+AL*AG 

200 continue 

201 IHsIDIAaClC) 

AK(XN)*AX(IN)-TR 

IF (AK(IN) .LT.ZERO) IDIAG<IC)— IDXAG(IC) 

ioo co st un i s 
RETURN 

901 WRITE (NV.lOOl) 

STOP 

1001 FORMAT (SX,29H«**ST0P ERROR IN ID I AG VECTOR) 

1010 FORMAT (SI, ’ CONDITIONING OF THE STIFFNESS MATRIX '/2X, 

1 'KIN D1AG TERMs' .1PE10.3,' XAZ DIAG TERH=' .E10.3) 

END 

cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc 

SUBROUTINE NULL (AK,AH s VL,VR,IDIAG,ISL,NDOF,NRHOD) 

C 

C PURPOSE : CALCULATE THE NULL SPACE OF AX AND PUT AN 
C ORTHONORMALXSES BASE OF THIS SPACE INTO THE 

C NRMQD FIRST VECTORS OF VR. 

C VL ■ AM x VR AFTER EXECUTION 

C 

C 

IMPLICIT REAL*8 (A-H.Q-Z) 

DIMENSION AK(1) ,AM(1) ,IDIAG(1) ,ISL(1) ,VL(NDOF,*) ,VR(NDOF,*) 

C 

C STORE THE SINGULAR COLUMNS INTO THE BEGINNING OF VR 
C 

C THE SINGULAR EQUATION ARE NOV INDEXED BY NEGATIVE VALUES 
C INTO ISL INSTEAD OF ID I AG 
C 

NRMQD >0 

DO 1 IC*1,ND0F 

IF (IDIAG(IC).GT.O) GOTO 1 
IDIAG(IC)a-IDXAG(XC) 

NRHOD*NRMQD+ 1 
MIC*ISL(IC) 

IH=IDIAG(IC)-IC 
DO 2 K»i,MIC-l 

2 VR(K,NRM0D)=0. 

DO 3 K-MIC.IC-l 

VR(X,NRMOD)»AX(IN+K) 

AK(IN*K)»0. 

3 CONTINUE 
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ISL(IC)«-ISL(IC) 

Va(IC,SHBOD)»-l. 

AK(XDXAG(IC)) a l. 

DO 4 X«ICn,HDOF 
4 VR(K ,HRHQD)*0 . 

1 COHTXHUE 
C 

C BAXSUBSTITUTE 

e 

DO 200 Hal.HBXOD 
IIC-HDOF 

DO 200 XC*2,HD0F 
TRaVR(IIC.H) 

IClaIIC-1 

IE«ISL(IIC) 

m»IDIAO(IIC)-IIC 

IF ((IX.GT.IC1) . OE. (IK.LE.O)) GO TO 221 
DO 220 K-XK.ICl 

220 COHTXHOB 

221 IICaIIC-1 
200 COHTZHOE 

C 

C ORTHOGOHALISATIOH , ; ; 

C 

DO .10 H»1,HRH0D 
DO 20 Kal.H-l 
TRaO . 

DO 30 lal.HDOF 

30 TR»TR*VI.(I,X)*VRU 1 H) 

DO 40 lal.HDOF , 

40 VR(I ,H)aVR(I.,H)“TRaVR(I ,X) 

20 COHTIHUE 

CALL HULT(AH,VL(1 ,H) ,VR(1 ,H) .ISL.IDIAG , 1 ,NDOF) 

TRaO. 

DO 50 lal.HDOF 

SO TRaTRm<I,H)*VL(I-,H) 

TRal/SQRT(TR) 

DO 60 Xal.HDOF 

VR(I,H)aVR(I.H)*TR 
VL(I,H)»VL(I,H)*TR 
60 COHTIHUE 

10 COHTIHUE 

RETURH 

EHD 

ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc 

SUBROUTXHE ORTHO ( VL, VR, HDOF, NRMOD ,HVEC) 

C 

C PURPOSE : ORTHOGOHALXSE THE HSOB LAST COLUKHS OF VL (LAST 
C EVALUATED SOLOTIOH) WITH RESPECT TO THE HULL SPACE 

C OF A. FOR THE AH SCALAR PRODUCT 

C 
C 

IMPLICIT REAL*8 (A-H.O-Z) , 

XHTEGER HDOF.HRHOD.HVEC 

DIKEHSIOH VLCNDOF ,HVEC) , VR(HDOF ,HVEC) 

HSUBaHVEC-HRMOD 
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DO 1 Jal.HSUB 
DO 2 2*i,HRSO0 
S-O. 

DO 3 Kat.IDOF 

3 S»S*VL(K , I) *VL(K .NRHOD* J ) 

DO 4 L»l t HDOF 

4 7L(L,HHJI0D+J)aVI,(L,HIUC0D+J)-S*Va(L,I) 

a coHTuroE 

1 COSTXHUE 

KSTOBV 

END 


ccccccccccecccccccccceccccccccccccccccccccccceccccccccccccccccccccccc 


File: animout.f 


C-Module ANXMOUT 

OAutbor K. Alvin 

C=Date June 1990 

C^Bloci; Fortran 

C 

C subroutine AHIHOOT 

C 

C Purpose: 

C This subroutine produces an output file to be used to 

C visualize tbe simulation using MESH. 

C 

C . 

c 

subroutine AHMOOT(q I id,nnp .time ,out) 
real*8 q(l) .time ,v(3) 
integer out,id(6,<),uup,i,k 
vrite(out,'(gl5.S)’) time 
do 100 isl.nnp 
do 200 k«l,3 

if (id(k,i) .ne. 0) then 
v(k) » q(id(k,i)) 
else 

v(k) » 0. 
end if 

200 continue 

write (out, 1000) (v(k),k*l,3) 

100 continue 

1000 format(3glS.S) 
return 
end 


File: stiffrc.f 


C subroutine STIFFRC ( r .fact ,kq) 
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racoraiva aubroutiaa STIFFRC(v,f act ,kq) 
rad*8 T(l),fact,kq(l) 
includa 'sharad.iac' 

C ASSEMBLE EACH ELEMEHT MASS AMD STIFFNESS 

do 100 ii*l .adoaain 

CVOt CHCALL 

do 110 jjal.aald(ii) 
a a alau»(jj,ii) 
call ELEFRC(7,*act,kq,a) 

110 eoatiaaa 

100 eoatiaaa 

ratora 

•ad 

C aabroatiao ELEFRC(v,tact ,kq,n) 

racuraiva anbrontiaa ELEFRC(7,f act ,kq,a) 

raal*8 7(1) ,fact,kq(l) 
iatagar a 

iaduda ’ aharad.iac' 

C LOCAL VARIABLES 

para»atar(HAXSEQ°24) 
raal*8 ak(MAISEQ .MAISEQ) 
iatagar lm(HAXSEQ) t aaaq 

do 20 k»l,4 
jaiz(k.a) 

it ((atypa(a) .aq.l) ,aad.(k.gt .2)) j » 0 
do 30 ial,6 
kkafla(k-l) ♦ i 
it (j .aa. 0) thaa 
la(kk) » idCi.j) 
alaa 

la(kk) > 0 
•adit 

30 eoatiaaa 

20 eoatiaaa 

acaqsl2 

call L0A0SK(ak,a,aaaq) 

call ESTIFVMCak.la.asaq.T.kq.fact) 

ratura 

•ad 
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C subroutine LOADSK(ak,a,naeq) 

recursive subroutine LOADSKCsk.n.nseq) 

include ’ shared. inc’ 

6 

real*8 sk(nseq,l) 
integer a, nseq 

k«0 

do 10 jal t nseq 
do 20 i«i,j 

k«kn 

ek(i, j)aestifm(k,n) 
sk(j.i)»sk(i.j) 

20 continue 

10 continue 

return 

end 


File: estifvm.f 


C subroutine ESTIFVH(ak,lni,aaeq, v.kq.fact) 

recursive subroutine ESTIFVH ( ak , la , na aq , v , kq , f act ) 
C ARGOKEHYS 

real«8 sk(nseq.l) ,v(l) ,kq(l) ,iact 
integer ln(l),nseq 

do 20 j • 1, nseq 
k ■ l*(j) 

if (k . eq. 0) goto 20 
do 10 i = 1, nseq 
m a lffi(i) 

if (a .eq. 0) goto 10 
kq(m) a kq(m) ♦ sk(i, j)*v(k)*fact 
10 continue 

20 continue 

return 

end 

CaEnd Fortran 


File: renum.f 


CaDECK REHUH 

C=PURP0SE- REHUHBERS THE GRID P0ZHTS TO KI51HIZE PROFILE STORAGE 
C=A0TH0a V S BELYZH and DOC HGU7EH 7-5-90 
C 

subroutine REHTJH 
C 

include ' shared. inc’ 
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Initial!*# vector 
maxtry*2*nap/ 3 ♦! 


ater»s»nnp*aaxtry 
do 22 j*i,nteras 
22 iadjcy(j)«0 

do 10 i*l,nap 
10 icoaat(i)«0 

do 1 i*l,nel 
nodea»ix(l,i) 
aodeb>ix(2,i) 

if ((aodoa .eq. 0).or.(aodeb . sq. 0)) go to 1 
ieount (nodea) «icouat (aodoa) *1 
ieount (nodab) *icouat (nodab ) > 1 
lamicouatfaodoa) 
ib*icount (nodab) 

if (ia.gt .maxtry .or. ib.gt. maxtry) go to 34S 
locate a (aodea-l)*maxtry+ia 
ladjcy(locate)*nodeb 
loeatea(nodab-l)*maxtry*ib 
iad jcy (locate) *nodaa 
1 contiaua 


ii a 0 

do 37 i*l,ntarma 

if ( iadjcy(i) .eq. 0 ) go to 37 
ii«ii+l 

iadjcyCii)aiadjcy(i) 

37 coatiaaa 


last*0 

do 2 i*l,nnp 

laat =las t ♦icoont (i ) 

2 coatiaaa 

jjoicooat(l) 
icouat(l)«l 
do 3 i*2 ,nnp*l 
kk«icount(i) 
i cooat ( i ) aicooat ( i- 1 ) ♦ j j 
jj»kk 

3 coatiaaa 

go to 666 
346 vrite(6,&55) 

566 format(2x, • arror in dimansion for KAZTE7 !! ') 
556 coatiaaa 
c ««*****«******«**a«*aaa 

call GZKRClKnnp.icoant.iadjcy.para.mask.xls) 

ratoxa 

and 

cxxmxxxxxxxxxmxx 


subroutine ganrca(naqns , xad j , ad j acy , porm .mask , zls ) 

c rafaraaca: computer solution of large sparse positive definite 

c systems, alaa gaorga ft josepb w-h lin 

c (prentive-hall.iac. .eaglesood cliffs, HI 07632) 

integer adjncy(l) ,aask(l) ,parm(l) ,zls(l) 
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iateger xadj (1) , ccsize, i.neqna ,nlvl,aus,root 
do 100 ial,neqns 

100 coatiaao 
nasal 

do 200 1*1 ,neqne 

i£( aaak(i).eq.0 ) go to 200 
rootsi 

call fnroot(root , xadj , ad jncy, mask, nlvl.xls ,perm(num) ) 
call rca(root, xadj .adjacy .mask, perm(aum) , ccaize.xls) 
aasenuBtccsize 
if(au».gt .neqae) go to 987 


200 coatiaao 

c pera(nev node) sold aode 

c now, aaakCold node)* nee aode 


987 coatiaao 

do 11 aew»l, neqae 
iold*pexs(aev) 

®&sk(iold)*ne* 

e arite(6,*) ’ iold,mask(iold) * ' ,iold,mask(iold) 

11 coatiaao 
retara 
ead 

cxxxxx.,zm%xxxxmxxxxxmxxxxmxxxxxxxxxxxxxxxxxxxxx 

subroutine faroot (root , xadj .adjncy, mask, nlvl.xls, Is) 
integer adjacy(l) ,la(l) ,mask(l) ,xls(l) 
iateger xadj(l) , ccsize, j , jatrt,k,kstop,kstrt, 

$ misdeg.nabor.ndeg.nlvl, node, nualvl, root 

call rootla(root,xadj , ad j acy , mask , alvl , xla ,1s) 
ceaizesxla(alvl+l)-»l 

if(nlvl.eq.l .or. alvl.eq. ccsize) return 
100 jatrt»xla(nlvl) 
miadegaccaixe 
root»la(j8trt) 

if ( ccsize. eq, jatrt) go to 400 
do 300 jsjstrt, ccsize 
nod*«Xs(j) 
adegaO 

kstrt»xad j (node) 
katopaxadj (node+l)-l 
do 200 kakatrt .kstop 
aabor=adjncy(k) 

if( mask(aabor) .gt. 0) ndeg=ndeg+l 
200 coatiaao 

if (adeg.ge.aindeg) go to 300 
root*node 
nindeg*ndeg 
300 coatiaao 

400 call rootla (root, xadj , ad j acy, mask, nualvl.xls ,1s) 
if (aaalvl.le.alvl) return 
alvl»nanlvl 

if (alvl.lt. ccaize) go to 100 

retara 

ead 

cxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxx 

aabroatiae rca(root,xadj .adjncy, mask, perm, ecaize.deg) 
iateger adjncy(l) ,deg(l) ,aask(l) ,perm(l) 
integer xadj(l) ,ccaize,fabr,i,j ,jstop,jstrt,k,l,lbegin, 
8 lnbr.lperm.lvlead.abr, node, root 
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call degr**(root,x&dj .adjney, mask, d*g,ecaixe, pens) 

aa*k(root)«0 

if (ccsixe.le.l) return 

lTl*nd*0 

lnbr*l 

100 lbagiaalvlandal 

lTl«ad«bbr 

do 600 i<*lb*gin,lvl*ad 
node*perm(i) 
j strt*xad j (nod*) 
j»top*xadj (nod*+l)-l 
fnbr*labr*l 
do 200 j*jstrt , jatop 
abr»adjney(j) 

if (maak(nbr) . eq.0) go to 200 

labr«labr*l 

Bask(nbr)>0 

p«ra( lnbr) ®nbr 

200 continue 

if(fnbr .g*.lnbr) go to 600 
k»£nbr 

300 l»k 

k»k*l 

nbr»pora(k) 

400 if(l.lt.fnbr) go to S00 

lp*ra*p*rm(l) 

if( d«g(lp«rm) . 1* . dog(nbr) ) go to S00 
p*rm(l*l)*lp*ra 

lBl_i , 
go to 400 

SOO p*na(l*l)°nbr 

if (k. It. lnbr) go to 300 

600 coot inn* 

if (lnbr.gt .lvlend) go to 100 
k*ccsiz*/2 
l*ccsixe 
do 700 i*l,k 
lp*rw*p*ra(l) 

P*na( 1) *p*r»( i) 

p*xm(i)*lp*ra 

1 » 1-1 

700 continue 

r*tnm 
•nd 

subroutine rootls(root ,xsdj , ad jncy.aaak, nlvl.xls ,1s) 
int*g*r adjncy(l) ,ls(l) ,mask(l) ,xls(l) 
integer xadj(l) ,i,j , jstop, jetrt.lbegia.ccsize.lvlend, 
$ lxsine.nbr ,alvl, node, root 

mask ( root )<*0 
la(l)*root 

nlrl»0 

lxl*nd«0 

eesiXMl 

200 lbegin*lvl*nd>l 
lrl«nd>eceize 
nlvlenlul+i 
xls ( nlvl) »lb*gin 
do 400 i*lb*gin, lvlend 
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aedasls(i) 
j strtax&dj (node) 
j stop*xad j (aod«> i ) -i 
12 (jatop.lt. jatrt) go to 400 
do 300 j* jatrt, jatop 
abroad jncy(j) 

if ( aaak(nbr) .aq.O) go to 300 
eeaizaseeaiza+l 
la(ceaiza)aabr 
B&sk(nbr)*0 
300 eoatiaaa 
400 eoatiaaa 

lvsiz**ccsiza-lvlaad 
if (lvaiza. gt .0) go to 200 
xls(nlvl*l)alvlead+l 
do 800 i»l,ecsiza 
noda«la(i) 
aaak(noda)*l 
500 Croatian* 

ratara 
and 

c mxzmmmama 

subroutisa dagr aa (root, xad j .adjncy.naak.deg.ccsiza.ls) 
intagar adjncy(l) ,dag(l) ,la(l) .oask(l) 
iatagar xadj(l) ,ccaiza,i,idag, j , jatop, jstrt, 

$ Ibagin , lvland , lvaiza . abr , node . root 

la(l)«root 

xad j (root ) "-xad j (root ) 

lvlaadaO 

ccaiza»l 

100 lbagia»lvland>l 
lvlaadaceaiza 
do 400 islbagia, lvland 
aoda»la(i) 
j a t rt =-xad j ( noda ) 
j»top*iab»( xadj(noda+l) ) -1 
idagaO 

if (jatop. It. jatrt) go to 300 
do 200 j® jatrt, jatop 
nbr®adjney(j) 

if( naak(abr) .aq.O ) go to 200 
idag*idag+l 

if (xad j (abr) .It. 0) go to 200 
xad j (nbr)a-xadj (abr) 
ccaiza*ccaiza+l 
la(ccaiza)®nbr 
200 eoatiaaa 

300 dag(noda)=idag 
400 eoatiaaa " 

lvaiza*ceaiza- lvland 
if (lvaiza. gt.0) go to 100 
do 500 i*l,ccaize 
noda®la(i) 

xadj (noda)s-zadj (node) 

500 eoatiaaa 

retarn 
and 


File: kfilter.f 
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C subroutine KFILTER 

recursive subroutine KFILTER 
include ’shared.inc* 

f 

do 10 i » l.ndof 

go(i) * delsq*(f (i)+go(i))+delta*pe(i)4aass(jdiag(i))*qe(i) 
gk(i) « delta*(f (i)*gk(i))+pe(i) 

10 continue 

call SOLYER(eo,go,jdiag,ndof ,2) 

Activate EBE computations for internal force by using STIFF8.C 
subroutine. Otherwise use PKVMAD (profile matriz/vector mult-add 

C call ?ff¥KAD(stii .jdiag, go, ndaf, -delta, gk,A.dO) 

call STIFFRC (go, -delta ,gk) 

do 100 i a l.ndof 

qe(i) o 2.*go(i) - qe(i) ' 

pe(i) » 2.*gk(i) - pe(i) 

100 continue 

return 

end 


<3 
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Abstract 

A general form for the first-order representation of the continuous, second-order linear 
structural dynamics equations is introduced in order to derive a corresponding form of 
first-order continuous Kalman filtering equations. Time integration of the resulting first- 
order Kalman filtering equations is carried out via a set of linear multistep integration 
formulas. It is shown that a judicious combined selection of computational paths and the 
undetermined matrices introduced in the general form of the first-order linear structural 
systems leads to a class of second-order discrete Kalman filtering equations involving only 
symmetric, sparse N x N solution matrices. The present integration procedure thus over- 
comes the difficulty in resolving the difference between the time derivative of the estimated 
displacement vector ( Jjx) and the estimated velocity vector (x) that axe encountered when 
one attempts first to eliminate (x) in order to form an equivalent set of second-order fil- 
tering equations in terms of (-^x). A partitioned solution procedure is then employed to 
exploit matrix symmetry and sparsity of the original second-order structural systems, thus 
realizing substantial computational simplicity heretofore thought difficult to achieve. 

+ An earlier version of the present paper without numerical experiments was presented at 
the AIAA Guidance and Control Conference, Port aland, Ore., 20-22 August 1990, Paper 
No. AIAA 90-3387. 
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Introduction 

Current practice in the design, modeling and analysis of flexible large space structures 
is by and large based on the finite element method and the associated software. The 
resulting discrete equations of motion for structures, both in terms of physical coordi- 
nates and of modal coordinates, are expressed in a second-order form. As a result, the 
structural engineering co mmuni ty has been investing a considerable amount of research 
and development resources to develop computer-oriented discrete modeling tools, analysis 
methods and interface capabilities with design synthesis procedures; all of these exploiting 
the characteristics, of second-order models. 

On the other hand, modern linear control theory has its roots firmly in a first-order form 
of the governing differential equations, e.g., (Kwakemaak and Sivan, 1972). Thus, several 
investigators have addressed the issues of interfacing second-order structural systems and 
control theory based on the first-order form (Hughes and Skelton, 1980; Arnold and Laub, 
1984; Bender and Laub, 1985; Oshman, Inman and Laub, 1987; Belvin and Park, 1989,, 
1990). As a result of these studies, it has become straightforward for one to synthesize 
non-observer based control laws within the framework of a first-order control theory and 
then to recast the resulting control laws in terms of the second-order structural systems. 

Unfortunately, controllers based on a first-order observer are difficult to express in a pure 
second-order form because the first-order observer implicitly incorporates an additional 
filter equation (Belvin and Park, 1989). However a recent work ( Juang and Maghami, 1990) 
has enabled the first-order observer gain matrices to be synthesized using only second-order 
equations. To complement the second-order gain synthesis, the objective of the present 
paper is to develop a second-order based simulation procedure for first-order obsen 3rs. 
The particular class of first-order observers chosen for study are the Kalman Filter based 
state estimators as applied to second-order structural systems. The procedure permits 
simulation of first-order observers with nearly the same solution procedure used for treating 
the structural dynamics equation. Hence, the reduced size of system matrices and the 
computational techniques that are tailored to sparse second-order structural systems may 
be employed. As will be shown, the procedure hinges on discrete time integration formulas 
to effectively reduce the continuous time Kalman Filter to a set of second-order difference 
equations. 

The paper first reviews of the conventional first-order representation of the continuous 
second-order structural equations of motion. An examination of the corresponding first- 
order Kalman filtering equations indicates that, due to the difference in the derivative of 
the estimated displacement (J^x) and the estimated velocity (i), transformation of the 
first-order observer into an equivalent second-order observer requires the time derivative 
of measurement data, a process not recommended for practical implementation. 
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Next, a transformation via a generalized momentum is introduced to recast the structural 
equations of motion in a general first-order setting. It is shown that discrete time numerical 
integration followed by reduction of the resulting difference equations circumvents the need 
for the time derivative of measurements to solve Kalman filtering equations in a second- 
order framework. Hence, the Kalman filter equations can be solved using a second-order 
solution software package. 

Subsequently, computer implementation aspects of the present second-order observer are 
presented. Several computational paths Eire discussed in the context of discrete said con- 
tinuous time simulation. For continuous time simulation, an equation augmentation is 
introduced to exploit the symmetry and sparcity of the attendant matrices by maintain- 
ing state dependant control Eind observer terms on the right-hand-side (RHS) of the filter 
equations. In addition, the computational efficiency of the present second order observer 
as compared to the first order observer is presented. 

Continuous Formulation of Observers 
for Structural Systems 

Linear, second-order discrete structural models can be expressed as 

Mx + Dx + Kx — Bu + Gw , x(0) = x 0 , x(0) = x 0 (1) 

u = —Z\X — Z 2 x 

with the associated measurements 

z = H\x + H 2 x + v * (2) 

where M,D,K are the mass damping and stiffness matrices of size (N x N)\ x is the 
structural displacement vector, ( N x 1); u is the active control force (m x 1); B is a 
constant force distribution matrix ( N x m); z is a set of measurements (r x 1); Hi and H 2 
are the measurement distribution matrices (r x iV); Z\ and Z 2 are the control feedback gain 
matrices (m x N)] w and v Eire zero-meem, white Gaussian processes with their respective 
covariances Q and R) and the superscript dot designates time differentiation. In the present 
study, we will restrict ourselves to the case wherein Q and R axe uncorrelated with each 
other and the initial conditions x 0 and x 0 are also themselves jointly Gaussian with known 
means and covariances. 

The conventional representation of (1) in a first-order form is facilitated by 

{ X\ = x 

x 2 = x = ii (3) 

Mx 2 = Mx = Bu + Gw — Dx 2 — Kx\ 
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which, when cast in a first-order form, can be expressed as 

{ 


Eq = Fq + Bu + Gw, q = ( x x x 2 ) T 
z — Hq + v 


( 4 ) 


where 


Mo y ■ m-°* ~ d \ • 

MM M£} 


(5) 


It is well-known that the K a.lma.rt filtering equations (Kalman, 1961; Kalman and Bucy, 
1963) for (4) can be shown to be (Arnold and Laub, 1984): 


Eq = Fq +.Bu + EPH T R~ l z 


where 


z — z — Hq , P = 


u S T 
S L 


A J Xj 1 J X 1 

1 X 2 J 1 X j 


( 6 ) 


(7) 


in which JJ and L are positive definite matrices and the matrix P is determined by the 
Riccati equation (Kwakemaak and Sivan, 1972; Arnold and Laub, 1984) 


EPE t = FPE t + EPF t - EPH t R~'HPE t + GQG t 


( 8 ) 


The inherent difficulty of reducing the first-order Kalman filtering equations given by (6) 
to second order form can be appreciated if one attempts to write (6) in a form introduced 
in (3): 

a) i\ = x 

b) x 2 = x = £\ — L\z (9) 

c) Mx 2 = —Dx 2 — Kx i + Bu + ML 2 z 


\ 


where 


L x = (H X U + H 2 S) t R-\ L 2 = (H X S T + H 2 Lf R 


Td-1 


Note from (9b) that x 2 i\. In other words, the time derivative of the estimated dis- 
placement (z) is not the same as the estimated velocity (z); hence, x x and z 2 must be 
treated as two independent variables, an important observation somehow overlooked in 
Hashernipour and Laub (1988). 

Of course, although not practical, one can eliminate x 2 from (9). Assuming x x and x 2 are 
differentiable, differentiate (9b) and multiply both sides by M to obtain 


Mx i = Mx 2 + ML\'z 


( 10 ) 
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Substituting Mi 2 from (9c) and X 2 from (9b) in (10) yields 

Mx\ = —D(x\ — L\z) — Kx\ + Bu + ML 2 Z + ML\i (11) 

which, upon rearrangements, becomes 

Mxi + Dxi + Kx\ = Bu + ML 22 + ML\z + DL\z (12) 

There are two difficulties with the above second-order observer. First, the numerical 
solution of (12) involves the computation of xi when rate measurements are made. The 
accuracy of this computation is in general very susceptible to errors caused in numerical 
differentiation of xi. Second, and most important, the numerical evaluation of z that 
is required in (12) assumes that the derivative of measurement information is available 
which should be avoided in practice. We now present a computational procedure that 
circumvents the need for computing measurement derivatives and that enables one to 
construct observers based on the second-order models. 


Second-Order Transformation of 
Continuous Kalman Filtering Equations 

This section presents a transformation of the continuous time first-order Kalman filter to 
a discrete time set of second-order difference equations for digital implementation. The 
procedure avoids the need for measurement derivative information. In addition, the spar- 
sity and symmetry of the original mass, damping and stiffness matrices can be maintained. 
Prior to describing the numerical integration procedure, a transformation based on gener- 
alized momenta is presented which is later used to improve computational efficiency of the 
equation solution. 


Generalized Momenta 


Instead of the conventional transformation (3) of the second-order structural system (1) 
into a first-order form, let us consider the following generalized momenta (Jensen, 1974; 
Felippa and Park, 1978): 


{ a) x\ — x 

b ) X2 = AMii+Cxi 


(13) 


where A and C are constant matrices to be determined. Time differentiation of (13b) 
yields 

X 2 = AMx 1 + Cx 1 (14) 
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Substituting (1) via (13a) into (14), one obtains 


x 2 — A(Bu + Gw) — (AD — C)x\ — AKx\ 


(15) 


Finally, pairing of (13b) and (15) gives the following first-order form: 

AM OlJx,! f C -I] fzil 
AD-C l\ \i 2 J [AK 0j\x 2 j 

[a(Hu + Gu>)] ^ 

The assoriated Kalman filtering equation can be shown to be of the following form: 

!][{;]■ <") 


AM 

AD-C 


where 


L x =(HxU + H 2 S) t RT 1 , L 2 = (H\S T + H 2 L) T R 


-(U.qT 


and H\ and H 2 correspond to a modified form of measurements expressed as 


z = H\X H 2 x = HxX\ + H 2 x 2 


(18) 


where 

Hx=Hx- H 2 M~ l A~ x C , H 2 = H 2 M~ l A~ l 

Clearly, as in the conventional first-order form (9), Xi and x 2 in (17) axe now two inde- 
pendent variables. Specifically, the case of A = M -1 and C = 0 corresponds to (3) with 
x 2 = xi. However, as we shall see below, the Kalman filtering equations based on the 
generalized momenta (13) offer several computational advantages over (3). 


Numerical Integration 

At this juncture it is noted that in the previous section one first performs the elimination 
of i\ in order to obtain a second-order observer, then performs the numerical solution 
of the resulting second-order observer. This approach has the disadvantage of having to 
deal with the time derivative of measurement data. To avoid this, we will first integrate 
numerically the associated Kalman filtering equation (17). 
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The direct time integration formula we propose to employ is a mid-point version of the 
trapezoidal rule: 



(19) 


where the superscript n denotes the discrete time interval t n = nh , h is the time increment 
and 6 = h/2. 

Time discretization of (17) by (19a) at the n + 1/2 time step yields 



The above difference equations require the solution of matrix equations of 2 N variables, 
namely, in terms of the two variables x £ +1 ^ 2 and x” +1/f2 , each with a size of N. To 
reduce the above coupled equations of order 2JV into the corresponding ones of order N, 
we proceed in the following way by exploiting the nature of parametric matrices of A and 
C as introduced in (13). To this end, we write out (20) as two coupled difference equations 
as follows: 

AM(i" +1/J - x”) + <5(C r Xj +1/2 - x£ +1/2 ) 

= 6AML x z n + 1 ! 2 (21) 

(AD - C)(x" +1/2 - x?) + (x£ +1/2 - x J) + 6AKx” +1/2 

= 6(AD - C)l x z n+l ! 2 + 6L 2 z n+l12 + 6ABu n+1 ' 2 (22) 

Multiplying (22) by 8 and adding the resulting equation to (21) yields 

A(M + 6D + 6 2 K)x”+ 1/2 = (AM + 6(AD - C))x ? + 8x^ 

+{8AML X + 6 2 (AD - C)L X + 6 2 L 2 }z n+1/2 + 6 2 ABu n+1 ' 2 (23) 


Of several possible choices for matrices A and B , we will examine 

f a) A = I , C = D 
\ b) A = M~ l , C = 0 


( 24 ) 
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The choice of (24a) reduces (23) to: 

(M + 8D + 8 2 K)x" +1/2 = Mi” + $x£ + 8 2 Bu n+1 ' 2 

+6{MLi+SL 2 }z n+1/2 (25) 

so that once i” +1 ^ 2 is computed, Xj"*" 1 / 2 i g obtained from (22) rewritten as 

x% +1/2 = x% + 6g n - 8Kx” +1/2 (26) 

where 

g n = Bu n+1 ^+L 2 z n+1/2 (27) 

which is already computed in order to construct the right-hand side of (25). Hence, 
_Kx™ +1 / 2 is the only additional computation needed to obtain x” +1 ^ 2 . It is noted that 
neither any numerical differentiation nor matrix inversion is required in computing x£ +1 ^ 2 . 
This has been achieved through the introduction of the general transformation (13) and 
the particular choice of the parameter matrices given by (24a). 

On the other hand, if one chooses the conventional representation (24b), the solution of 
-n+i /2 - s Stained from (23) 


(M + 6D + S 2 K )x x f 1/2 = (M + 6D)x ” + 8Mx^ 

+8{(M + 8D)L 1 +8ML 2 }z n+1/2 + 6 2 Bu n + 1 ' 2 (28) 

Once x” +1 / 2 is obtained, x£ +1 ^ 2 can be computed either by 

x” +1/2 = (x" +1/2 - x”)/8 - L x z n+1 / 2 (29) 


which is not accurate due to the numerical differentiation to obtain x x + ^ , or by (22) 


x” +1/2 = xj + 8g n - 8M~ l Kx” +l/2 - 

M~ 1 D(x” +1/2 — x”) + 8M~ 1 DLiz n+1 f 2 (30) 

which involves two additional matrix-vector multiplications, when D ^ 0, as compared 
with the choice of A = I and C = D . . Thus (24a) is the preferred representation in a 
first-order form of the second-order structural dynamics equations (1) and is used in the 
remainder of this work. 


8 


Decoupling Of Difference Equations 

We have seen in the previous section, instead of solving the first-order Kalman filtering 
equations of 2 n variables for the structural dynamics systems (1), the solution of the im- 
plicit time-discrete observer equation (25) of n variables can potentially offer a substantial 
computational saving by exploiting the reduced size and sparsity of M , D and K. This 
assumes that f n+1 / 2 and u n+1 / 2 are available, which is not the case since at the n th time 
step 


u n+l/2 = _Z X x n + 112 - Z 2 x” +1/2 

(31) 

* 

f n+l/2 = z „+l/2 _ tf^+1/2 _ £ 2f "+l/2 

(32) 


requires both x" +1 ^ 2 and x£ +1/ ^ 2 even if x n+1 / 2 is assumed to be known from measurements 
or by solution of (1). Note in (32), the control gain matrices are transformed by 

Z l = Z l - Z 2 M~ 1 A~ 1 C , Z 2 = Z 2 M~ x A~ l 

-*r 

There are two distinct approaches to uncouple (25) and (26) as described in the following 
sections. 

Discrete Time Update 

Equations (31) and (32) can be approximated using 

r +1/2 ~ z n _ ^ f n _ g 2i n (33) 

u "+i/2 ~-Z x x\-Z 2 x* (34) 

This approximation leads to a discrete time update of the control force and state correction 
terms which is analogous to that which exists in experiments where a finite bandwidth of 
measurement updates occurs. For discrete time approximation, the step size h = < n+1 — t n 
should be chosen to match the time required to acquire, process and output a control 
update. 

Discrete time simulation is quite simple to implement as the control force and state cor- 
rections are treated with no approximation on the right-hand-side (RHS) of (25) and (26). 
Should continuous time simulation be required, a different approach is necessary. 

Continuous Time Update 

To simulate the system given in (25) and (26) in continuous time, strictly speaking, one 
must rearrange (25) and (26) so that the terms involving x" +1/ ^ 2 and x£ +1 ^ 2 are augmented 
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to the left-hand-side (LHS) of the equations. However, this augmentation into the solution 
matrix ( M +6D+6 2 K ) would destroy the computational advantages of the matrix sparcity 
and symmetry. Thus, a partitioned solution procedure has been developed for continuous 
time simulation as described in (Park and Belvin, 1991). The procedure, briefly outlined 
herein, maintains the control force and state correction on the RHS of the equations as 
follows. 

First, x? 1 ' 2 and x£ +1/ ^ 2 are predicted by 


*"+1/2 _ 

z i p 


= x 


i > 


*"+1/2 _ -r 

X 2p ~ Z 2 


(35) 


However, instead of direct substitution of the above predicted quantity to obtain tt£ +1 ^ 2 
and 2 p +1 ^ 2 based on (31) and (32), equation augmentations are introduced to improve the 
accuracy of Up"*’ 1 ^ 2 and 2 p +1 ^ 2 . Of several augmentation procedures that axe applicable 
to construct discrete filters for the computations of u n + 1 / 2 and z"+i/ 2 , we substitute (26) 
into (31) and (32) to obtain 


u n+i/2 _ -z x x*+ 1/2 - Z 2 (x” - 8Kx" +1/2 + 
6Bu n+1 ' 2 +6L 2 z n+1 ' 2 ) 
r+l/2 »*«+!/* -,ffi*?+l/2 — 

tf 2 (xj - 8Kx* +1/ 2 + 8Bu n + »/? + 6L 2 z n+1/2 ) 


(36) 


Rearranging the above coupled equations, one obtains 


'( I + 6Z 2 B ) 6Z 2 L 2 


6H 2 B ( I + SH 2 L 2 ) ] { *" +1 / 2 / 

-Z 2 x? - (Zi - 8Z 2 K)x" +1/2 


u n +l/2 ^ _ 




+1 /2 - H 2 x 2 - (Hi - 5H 2 AT)x” +1/2 


} 


(37) 


which corresponds to a first order filter to reduce the errors in computing x 2 = Mi + Dx. 
A second-order discrete filter for computing u and z can be obtained by differentiating u 
and z to obtain 

li = —ZiXi — Z 2 x 2 
z = i — Hiii — H 2 x 2 


{s 


(38) 


and then substituting xi and x 2 from (17). Subsequently, (19) is applied to integrate the 
equations for u and z which yields 

\I+6Z 2 B + 8 2 Z 1 M- 1 B 6(Z 2 L 2 + ZiLi-¥8ZiM- 1 L 2 ) 1 / u n+1/2 1 __ 

[ 8(H 2 B + 8HiM~ x B) I + 8Hi(Li + 8M~ l L 2 ) + 8H 2 L 2 J \ z" +1 / 2 J “ 

/«*!_«/ Zi M-'(*; - SKx" +,/2 - Dx"* 1 ' 1 ) + ZiKx" + ' /2 I f 0 ) 

tw 1 HiM- x (xZ -8Kx" +ll2 Dx” +1/2 ) + H 2 Kx” +1/2 J + \ z n+1/2 - z” J } 
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The net effects of this augmentation are to filter out the errors committed in estimating 
both x\ and x 2 . Solution of (39) for ti n+1 / 2 and f n+1 / 2 permits (25) and (26) to be solved 
in continuous time for x" +1 ^ 2 and x”" 1 " 1 ^ 2 . Subsequently, (29b) is used for x" +1 and x£ +1 . 

The preceding augmentation (39) leads to an accurate estimate of the control force and 
observer error correction at the (n+1/2) time step. Although (39) involves the solution of 
an additional algebraic equation, the equation size is relatively small ( size = number of 
actuators (m) plus the number of measurements (r) ). Thus, (39) is an efficient method 
for continuous time simulation of the Kalman filter equations provided the size of (39) is 
significantly lower than the first order form of (4). The next section discusses the relative 
efficiency of the present method and the conventional first order solution. More details on 
the equation augmentation procedure (39) may be found in Park and Belvin (1991). 

Finally, it is noted that by following a similar time discretization procedure adopted for 
computing x" +1 ^ 2 and x""*” 1 ^ 2 , the structural dyn ami cs equation (1) can be solved by 

/ (M + 8D + 8 2 K)x^ l/2 = Mi" + 8x% + 8 2 Bu n+1 ' 2 . . . 

\ X 2 +1/2 = x£ + 8Bu n + x ! 2 - 8Kx"+ 1/2 1 * 

Thus, numerical solutions of the structural dynamics equation (1) and the observer equa- 
tion (20) can be carried out within the second-order solution context, thus realizing sub- 
stantial computational simplicity compared with the solution of first-order systems of equa- 
tions (4) and the corresponding first-order observer equations (6). 

It is emphasized that the solutions of both the structural displacement x and the re- 
constructed displacement x employ the same solution matrix, (M + 8D + 8 2 K). The 
computational stability of the present procedure can be examined as investigated in Park 
(1980) and Park and Felippa (1983, 1984). The result, when applied to the present case, 
can be stated as 

£ 2 Amax < 1 (41) 

where A m ax is the maximum eigenvalue of 

(\ 2 I + \Z 2 B + Z 1 M- 1 B)y = 0 (42) 

Experience has shown that |A max | is several orders of magnitude smaller than /i m ax of 

the structural dynamics eigenvalue problem: 

fiMy = Ky (43) 

Considering that a typical explicit algorithm has its stability limit Umax ■ h < 2, the 
maximum step size allowed by (42) is in fact several orders of magnitude larger than 
allowed by any explicit algorithm. 
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Computational Efficiency 

Solution of the Kalman filtering equations in second-order form is prompted by the po- 
tential gain in computational efficiency due to the beneficial nature of matrix sparcity and 
symmetry in the solution matrix of the second-order observer equations. There is an over- 
head to be paid for the present second-order procedure, that is, the additional computations 
introduced to minimize the control force and observer error terms on the right-hand-side of 
the resulting discrete equations. The following paragraphs show the second-order solution 
is most advantageous for observer models with sparse coefficient matrices M, D and K. 

Solution of the first order Kalman filter equation (6) or the second-order form (25-26, 39) 
may be performed using a time discretization as given by (19). For linear time invariant 
(LTI) systems, the solution matrix is decomposed once and subsequently upper and lower 
triangular system solutions are performed to compute the observer state at each time step. 
Thus, the computations required at each time step result from calculation of the RHS 
and subsequent triangular system solutions. For the results that follow, the number of 
floating point operations per second (flops) are estimated for LTI systems of order 0(N). 
In addition, it is assumed that the mass, damping and stiffness matrices ( M,D and K) 
are symmetric and banded with bandwidth aN, where 0>< a < (0.5 — 277). 

The first-order Kalman filter equation (6) requires (AN 2 + 2 Nr + O(N)) flops at each 
time step. The discrete time second-order Kalman filter solution (25-26, 33-34) require 
(8a 2 IV 2 -f 2alV 2 +3Nm+ANr + 0(N)) flops and the continuous time second-order Kalman 
filter (25-26, 39) require (8a N 2 + 2 aJV 2 + 5 Nm + 6 Nr + (r + m) 2 + O(N)) flops at each 
time step. To examine the relative efficiency of the first-order and second-order forms, 
several cases are presented as follows. 

First, a worst case condition is examined whereby M, D and K are fully populated (a = 
0.5 — 277) and r = m = N. For this condition, the number of flops are: 

First Order 61V 2 4- O(N) 

'< Second Order Discrete 10 N 2 + O(N) 

k Second Order Continuous 181V 2 4- O(N) 

Thus, for non-sparse systems with large numbers of sensors and actuators relative to the 
system order, the first order Kalman filter is 300 percent more efficient than the second- 
order continuous Kalman filter solution presented herein. 

For structural systems, M, and K are almost always banded. In addition, the number 
of sensors and actuators is usually small compared to the system order N. Hence, the 
value of a for which the second-order form becomes more computationally attractive than 
the first order form must be determined. If the assumption is made that the number of 
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actuators (m) and the number of measurements (r) is proportional to the bandwidth ( 
r = m = alV), the value of a which renders the second-order solution more efficient is 
readily obtained. For the 'second-order discrete Kalman filter, when a < 0.394 the second- 
order form is more efficient. Similarly, the second-order continuous K alman filter form is 
more efficient when a < 0.279. Since a obtains values approaching 0 when a modal based 
structural representation is used with few sensors and actuators, the second-order form 
can be substantially more efficient than the classical first-order form. A more detailed 
discussion can be found in Belvin (1989). 


Implementation and Numerical Evaluations 

The second-order discrete K alman filtering equation derived in (25) and (26) have been 
implemented along with the stabilized form of the controller u and the filtered measure- 
ments z in such a way the observer computational module can be interfaced with the 
partitioned control-structure interaction simulation package developed previously (Belvin, 
1989; Belvin Park, 1991; Alvin and Park, 1991). Table 1 contrasts the present CSI simula- 
tion procedure to conventional procedures. It is emphasized that the solution procedure of 
the present second-order discrete Kalman filtering equations (25) and (26) follows exactly 
the same steps as required in the solution of symmetric, sparse structural systems (or the 
plant dynamics in the jargon of control). It is this attribute that makes the present discrete 
observer attractive from the simulation viewpoint. 

The first example is a truss beam shown in Fig. 1, consisting of 8 bays with nodes 1 and 
2 fixed for cantilevered motions. The locations of actuator and sensor applications as well 
as their directions are given in Table 2. Figures 2, 3 and 4 are the vertical displacement 
histories at node 9 for open-loop, direct output feedback, and dynamically compensated 
feedback cases, respectively. Note the effectiveness of the dynamically compensated feed- 
back case by the present second-order discrete Kalman filtering equations as compared 
with the direct output feedback cases. Figure 5 illustrates a testbed evolutionary model 
of an Earth-pointing satellite. Eighteen actuators and 18 sensors are applied to the sys- 
tem for vibration control and their locations are provided in Tables 3 and 4. Figures 6, 
7, and 8 are a representative of the responses for open-loop, direct output feedback, and 
dynamically compensated cases, respectively. Note that u, response by the dynamically 
compensated case does drift away initially even though the settling time is about the same 
as that by the direct output feedback case. However, the sensor output axe assumed to 
be noise-free in these two numerical experiemnts. Although the objective of the present 
paper is to establish the computational effectiveness of the second-order discrete Kalman 
filtering equations, we conjecture that for noise-contaminated sensor output for which one 
would apply dynamic compensated strategies, the relative control performance may turn 
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out to be the opposite. Further simulations with the present procedure should shed light 
on the performance of dynamically compensated feedback systems for large-scale systems 
as they axe computationally more feasible than heretofore possible. 

Table 5 illustrates the computational overhead associated with the direct output feedback 
vs. the use of a dynamic compensation scheme by the output present Kalman filtering 
equations. In the numerical experiments reported herein, we have relied on Matlab software 
package (Wolfram, 1988) for the synthesis of both the control law gains and the discrete 
Kalman filter gain matrices. It is seen that the use of the present second-order discrete 
Kalman filtering equations for constructing dynamically compensated control laws adds 
computational overhead, only an equivalent of open-loop transient analysis of symmetric 
sparse systems of order N instead of 2 N x 2 N dense systems. 


Summary 

The present paper has addressed the advantageous features of employing the same direct 
time integration algorithm for solving the structural dynamics equations also to integrate 
the associated continuous Kalman filtering equations. The time discretization of the re- 
sulting Kalman filtering equations is further facilitated by employing a canonical first-order 
form via a generalized momenta. When used in conjunction with the previously developed 
stabilized form of control laws (Park and Belvin, 1991), the present procedure offers a sub- 
stantial computational advantage over the solution methods based on a first-order form 
when computing with large and sparse observer models. 

Computational stability of the present solution method for the observer equation has been 
assessed based on the stability analysis result of partitioned solution procedures (Park, 
1980). To obtain a sharper estimate of the stable step size, a more rigorous computational 
stability analysis is being carried out and will be reported in the future. 
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Structure: 


o) 


where 


and 


Sensor Output: b ) 

Estimator: c) 

Control Force: d) 

Estimation Error: e) 


Mq + Dq + Kq = f + Bu + Gw 

q(0) = qo, q(0) = qo 
z = Hx + v 

Mq + Dq + Kq = f + Bu + ML27 
q(o) = 0 , q(0) = 0 

u + F2M -1 Bu = F2(M -1 p + L27) + Fiq 
7 + H v L 2l = z - - Bu) - H d 4 


Table la Partitioned Control-Structure Interaction Equations 


Structure: 




Sensor Output: b) 

Estimator: c) 


Control Force: d ) 

l Estimation Error: e) 


x = Ax + Ef + Bu + Gw 
q(0) = q 0 j q(0) - q 0 
z = Hx + v 

i = Ax + Ef + Bu + L7 
x(0) = 0 
u = -Fx 

7 = z - (H<iq + H„q) 


H = [H<i H„], L 


-{?}■ ■ -m 
■si- ■ 


o' 

~ M _1 _ 

, F = [Fi F 2 ] 
Table lb Conventional Control-Structure Interaction Equatioons 


A = 


r 0 

1 




0 

1 

s 

1 

-M -1 D 


B = 

M_1 B 
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TABLE 2a: 

Actuator Placement for Truss Example Problem 


Actuator 

Node 

Component 

1 

2 

y 

2 

18 

y 

3 

9 

y 

. 4 

9 

X 


TABLE 2b: 

Sensor Placement for Truss Example Problem 


Sensor 

Type 

Node 

Component 

1 

Rate 

2 

y 

2 

Rate 

18 

y 

3 

Rate 

9 

y 

4 

Rate 

9 

X 

5 

Position 

9 

y 

6 

Position 

9 

X 
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TABLE 3: 

Actuator Placement for EPS Example Problem 


Actuator 

Node 

Component 

1 

97 

X 

2 

97 

z 

3 

96 

X 

4 

96 

z 

5 

65 

y 

6 

68 

y 

' 7 

59 

y 

8 

62 

y 

9 

45 

y . 

10 

45 

z 

11 

70 

y 

12 

70 

z 

13 

95 

X 

14 

95 

y 

15 

95 

z 

16 

95 

<f>T 

17 

95 

<t>y 

18 

95 
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TABLE 4: 

Sensor Placement for EPS Example Problem 


Sensor 

Type 

Node 

Component 

1 

Rate 

97 

X 

2 

Rate 

97 

z 

3 

Rate 

96 

X 

4 

Rate 

96 

z 

5 

Rate 

65 

y 

6 

Rate 

68 

y 

7 

Rate 

59 

y 

8 

Rate 

62 

y 

9 

Rate 

45 

y 

10 

Rate 

45 

z 

11 

Rate 

70 

y * 

12 

Rate 

70 

z 

13 

Position 

95 

X 

14 

Position 

95 

y 

15 

Position 

95 

z 

16 

Position 

95 

<t>x 

17 

Position 

95 

4>y 

18 

Position 

95 

<t>Z 
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TABLE 5: 

CPU Results for ACSIS Sequential and Parallel Versions 


Model 

Problem 

Type 

Sequential 

Parallel 

54 DOF 

Transient 

4.5 


Truss 

FSFB 

9.4 



K. Filter 

13.0 

10.7 

582 DOF 

Transient 

98.6 

100.3 

EPS7 

FSFB 

190.2 

294.5 


K. Filter 

284.2 

321.5 
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Truss Model: Open Loop Transient Response 

6.0 x 10“ 4 
4.8 x 10 -4 
3.6 x 10~ 4 
2.4 x 10~ 4 
1.2 X 10~ 4 


Time, sec 

Node 9, uy 

Figure 2: Truss Transient Response 
Truss Model: Full State Feedback Response 

6.0 x 10~ 4 
4.8 x 10“ 4 
3.6 x 10~ 4 
2.4 x 10~ 4 
1.2 X 10~ 4 
0.000 

0.000 0.200 0.400 0.600 0.800 1. 

Time, sec 




Node 9, uy 

Figure 3: Truss FSFB Response 




Deflection 


Truss Model: Controlled Response with Kalman Filter 
6.0 x 10“ 4 

4.8 x 10 -4 

3.6 x 10“ 4 

2.4 x 10~ 4 

1.2 x 10“ 4 

0.000 

0.000 0.200 0.400 0.600 0.800 1.000 

Time, sec 

Node 9, uy 

Figure 4: Truss Response with Filter 
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EPS7 Model: Open Loop Transient Response 

2.0 x 10” 4 

) 

1.4 x 10“ 4 

8.0 x 10 -5 

2.0 x 10 -5 

4.0 x 10“ 5 

1.0 x 10 -4 

0.000 2.000 4.000 6.000 8.000 10 

Time, sec 

— Node 45, ux Node 45, uz 

Node 45, uy 

Figure 6: EPS Transient Response 
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EPS7 Model: Full State Feedback Response 

2.0 x 10 -4 
1.4 x 10~ 4 

8.0 x 10“ 5 

2.0 x 10~ 5 
-4.0 x 10“ 5 
-1.0 x 10“ 4 

0.000 2.000 4.000 6.000 8.000 10.000 

Time, sec 

Node 45, ux Node 45, uz 

Node 45, uy 

Figure 7: EPS FSFB Response 
Model: Controlled Response w/Kalman Filter 

2.0 x 10 -4 
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-4.0 x 10“ 5 
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Figure 8: EPS Response with Filter 
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ABSTRACT 

The equations of motion for structures with adaptive elements for vibration control 
are presented for parallel computations to be used as a software package for read- time 
control of flexible space structures. A brief introduction of the state-of-tb vart parallel 
computational capability is also presented. Time marching strategies are developed for 
an effective use of massive parallel mapping, partitioning and the necessary arithmetic 
operations. An example is offered for the simulation of control-structure interaction on 
a parallel computer and the impact of the approach presented herein for applications 
in other disciplines than aerospace industry is assessed. 


1. Introduction 

Active suppression of structural vibrations or active control of flexible structures has 
made considerable progress in recent years. As a result, it is now possible to actively 
suppress vibrations in mechanical systems emanating from machine foundations, in 
robotic manufacturing arms, truss-space structures and automobile suspension sys- 
tems. A common characteristic to these applications of active control theory has been 
its discrete actuators and discrete sensors, ranging from proof mass actuators and gyro 
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dampers to strain gages and accelerometers. Because most available discrete actuators 
are inertia force-oriented devices, actuation often triggers coupling between the actu- 
ator dynamics and structural transients. A practical consequence of such coupling is 
a limitation of achievable final residual vibration level if both the actuator and struc- 
ture possess insufficient passive damping leveL It is noted that structures made of 
high stiffness composite materials have very low intrinsic damping, hence limiting the 
achievable residual vibration level for space maneuvering and space disturbance rejec- 
tion purposes. This has been a motivating factor for the development of distributed 
actuators and sensors which are often embedded as an integral part of the structure 
so that control force can be effectively maintained by strain actuation, thus alleviating 
the undesirable actuator dynamics associated with inertia-force actuation. 

Various activities that are being pursued by many investigators on the subject of 
adaptive structures may be categorized into three major thrusts: device developments, 
control laws synthesis and experimental demonstrations, and hardware/software im- 
plementation. The device developments effort has been the objective of many material 
scientists (1-3). As the applications needs increase it is expected that functionally more 
reliable electrostrictive and magnetostrictive elements will be available for use in active 
control/strain damping with improved product quality. 

The study of control laws synthesis and demonstration employing adaptive ele- 
ments has been one of the predominant activities in recent years. As scientists accu- 
mulate experience in the characterization of the coupling between the structure and 
the adaptive element, the applications will then be expanded from the current beam- 
like structures to the truss long beams, plates and shells. In order to effectively uti- 
lize as many adaptive elements as necessary for actively controlling the vibration of 
such large-scale structures in real-time operations, it will be imperative that the soft- 
ware/hardware components in the real-time control loop must be able to process data 
fast enough so that control commands and the measurements can be carried without 
saturating and/or jamming the control system. 

With the advent of new technology in distributed actuators and sensors [4-9], it ap- 
pears that a combination of decentralized/ distributed and hierarchical control strategies 
can be a viable alternative to conventional centralized control strategies. The real-time 
computer control of such systems as well as design of such control systems through 
iterating on simulations and hardware realizations thus will require the processing of a 
vast amount of data from and to the distributed actuators and sensors. A significant 
part of such data processing for the decentralized actuators and sensors is planned to be 
self-managed, viz., there will be embedded microprocessors for each actuator and sensor 
pair or for each group of them. However, the necessary links between the decentralized 
control systems and the global control system as well as the necessary global control 
strategy will still require computational power far in excess of presently available real- 
time data processing capability. In addition, if one contemplates the performance of 
neural-network control or adaptive control for onboard real-time control of large-scale 
space structures, the computational need will dramatically increase beyond the current 
capability. As a case in point, even for the control of 20-bay truss beam vibrations by 
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three proof mass actuators and six sensors, NASA/Langley is relying on CRAY-XMP 
for adequate real-time data processing requirements. 

The objective of this paper is thus to present a computational framework by which 
one can bring the two emerging new technologies together, namely, the distributed actu- 
ators and sensors and the parallel computing capability, toward the real-time control of 
vibrations in large structural systems such as space stations, space cranes and in-space 
construction facilities. We will then discuss the potential for applying such a space 
technology to mitigate and/or minimize the earthquake damage of ground structures 
such as high-rise buildings, bridges and lifeline equipment. 


2. Models for Structures with Embedded Actuators and Sensors 

The coupling between the structural behavior and an adaptive electrostrictive element, 
whether it is embedded or surface-mounted, is primarily due to the following constitu- 
tive relation [3,10-12]: 


where e and v are the electrical displacement (charges/unit area) and the electric field 
(volt /unit area), <r and e are the stress and strain, and s, g and c are the constitutive 
coefficient matrices, respectively. For magnetostrictive elements, one needs to replace 
e and v by the magnetic field (H) and the magnetic induction (B), respectively, and 
the subsequent derivations will hold without any loss of generality. 



The coupled equations of motion for the structure and the adaptive elements 
can proceed by augmenting the standard procedure for the structure with the elec- 
tric transient equations plus the appropriate modification of the structural equilibrium 
equations that reflect the coupled constitutive equations (1). The resulting coupled 
structural-piezoelectric equations of motion take the following form [13-15]: 


( Structure: 


a) Mq + Dq + (K, + K„)q = f + Sa 


q(o) = q 0 . q(0) = q 0 


Sensor Output: 

Actuator: 

Controller: 


6) y = H p q 4- H r q t H a a 

c) a+0a = B a u-S r |^ 

d) u + Gu = Ly 


(2) 




where 
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In the preceding equations, M is the mass matrix, D is the damping matrix, K, 
is the stiffness matrix due to structural strain-displacement relations and K a is the 
stiffness matrix due to the strain actuation. f(t) is the applied force. S is the actuator 
projection matrix. H p , H r and H, are the sensor calibration gain matrices, © is the 
actuator dynamic characteristics, B„ is the gain matrix that translates the applied 
current /charge and voltage into the corresponding strain and strain rate where S is 
the transducer conversion gain, q is the generalized displacement vector and and the 
superscript dot denotes time differentiation, and u is the control law that consists of 
the applied current (or charge), Io, and voltage across the electrostrictive devices, Vo, 
G is the electric circuit characteristics, and L is the optimum direct feedback gain 
matrix. The case of dynamic compensations can be augmented to (2) by introducing 
an observer. But in subsequent discussions we limit ourselves to direct feedback cases 
only. 


It is noted that the control laws, unlike conventional control-structure interaction 
systems, are not directly fed back into the structural equations. Instead, the controller 
is simply a regulator controlling the electric charge, the voltage or the current. These 
regulated electric quantities are then fed into the piezoelectric sensors and actuators. 
Hence, it is the piezoelectric actuation that triggers feedback into the structures. 


3. Parallel Computations for the Dynamics of Adaptive Structures 

The earliest recorded computational results in mechanics were the parabolic trajectory 
calculations of a falling body by Galileo [16]. Since then, most scientific computations 
have been carried out by anthropomorphic algorithms, viz., step-by-step binary and/or 
decimal arithmetics. To set the stage properly for the present objective, parallel com- 
putations of the dynamic response of structures with distributed adaptive elements, we 
recall a passage by Kepler ,to John Napier, the inventor of a logarithmic table: 

Newton was essentially dependent upon the results of Kepler’s cal- 
culations, and these calculations might not have been completed but 
for the aid of that logarithms afforded. Without the logarithms, ...-, . 
the development of modern science might have been very different 

( 17 ]. 

In terms of the present day data processing requirement, Napier’s logarithmic 
table in 1614 contained about 100 kilobytes, which was perhaps the most important 
computational aid to Kepler and Newton. Three and one-half centuries later we are 
witnessing gigabytes of tables being stored and retrieved at our disposal [18]. But 
these tables complement the weakness of the human mind and computational speed: 
long term memory and human arithmetic speed. In addition, for problems requiring a 
sequential nature of computations, i.e., ballistic trajectories which deal only with the 
position and velocity of a single shell or quasi-static equilibrium equations of a building 
structure, the computing activities do not interact with “time” and the computing 
efficiency affects only the humanpower efficiency for completing the computational 
task. . 
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There are many important scientific and engineering endeavors whose computa- 
tions must be fast enough for real-time delivery of the computed results. A classical 
example was Richardson's lattice model for weather prediction by numerical process 
in 1922 [19]. The motivation for adopting such a lattice concept was due to the fact 
that the equation state at each lattice node takes on a different value set in time and 
an efficient way of interchanging and transmitting the nodal values at each time step 
was mandatory if the computations were to be carried out in real-time to predict the 
weather. Indeed, this was the dawn of the parallel computing era, even though the basic 
idea had to wait for its validity for 60 years. Today, many controls engineering activi- 
ties have been implemented by using computers so that their intended functions can be 
monitored and controlled in real-time. These include chemical processing, autopiloting 
and vibration control of simple structures. It is important to note that the computa- 
tional framework employed for such applications is based on sequential architecture. 
Hence, we believe that future improvements that can deal with large parameter models 
and large parameter controls must adopt a parallel computational framework. One 
such area is the dynamics and control of large structures with distributed/embedded 
adaptive elements. 

In order to carry out the necessary parallel computations, there are three distinct 
steps that must be addressed: discretizing the structure into appropriate partitions, 
mapping the physical partitions onto the processors, and step advancing of the equation 
states. These will be discussed below. 

3.1 Partitioning and Mapping of Adaptive Structures 

Ideally, if the sensor and actuator leads fall on the discrete nodes, no spatial interpola- 
tion would be necessary. However, such a situation is either difficult to realize or may 
prohibit the use of spatially convolving sensors [20] that axe known to filter certain 
harmonic signals for minimizing phase lag in the feedback loop: Hence, we will assume 
that the sensor and actuator characteristics can be interpolated to the discrete nodes; 
in this way the partition boundaries can be chosen arbitrarily regardless of the physical 
locations of the sensor and actuator leads. In addition, this approach can lead to a 
natural embedding of the sensor and actuator characteristics into the finite element or 
boundary integral structural models. Once the partitioning is accomplished, the next 
step is to map the discrete partitions for adaptive elements onto the corresponding 
multiprocessors. 

Consider an adaptive structure that has been modeled as a set of discrete elements 
as shown in Fig. 1. In a sequential computing environment, in order to advance the 
necessary computations for the present states, the arithmetic operations are carried 
out step-by-step for each node at a time. Hence, each nodal-state computations is 
performed in a manner similar to one courier delivering and picking up all the mails 
throughout the entire routes. In a parallel computational environment, in contrast, 
there can be as many couriers as necessary who comb through the routes concurrently 
in order to pick up and deliver till the mail at once. One of the most popular concepts 
in executing such tasks is the hypercube architecture (see Fig. 2) whose every node 
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is associated with a processor. Thus, to process the necessary computations for an 
adaptive structure with 19 partitions, one can assign the 19 adaptive elements to 19 
processors as shown in Fig. 3. The procedure for assigning the physical domain (ele- 
ments) to the parallel processors with minimal interprocessor communications is called 
mapping. 

Of several techniques available for the processor mapping of the computational 
domains [22], we will adopt a heuristic mapper developed by Farhat [23] since it can 
accommodate both the synchronous and asynchronous cases with robust and accept- 
able complexities. An application, of this mapping technique for modeling a bulkhead 
substructure for massively parallel computing is shown in Fig. 4. A similar mapping 
can be used for parallel computations of adaptive structures. 

3.2 Parallel Data Structure and Algorithms 

We will assume that each processor is assigned to carry out all the necessary computa- 
tions for at least one set of a sensor, an actuator, and a controller or a group of them. 
Therefore, the word partition does not necessarily imply a finite element: it can be a 
substructure, an element or even a sublayer within the composite layer that includes a 
sensor or an actuator. In carrying out the step- advancing in time, one may invoke an 
implicit or explicit direct time integration algorithm. When an implicit algorithm is 
employed, one needs to communicate not only the state variable vectors but also the 
associated matrices, i.e., the stiffness matrix, among the processors. Although we will 
show our results using implicit algorithms, we will, for illustrative purposes, restrict 
ourselves to an explicit direct time integration algorithm as it is intrinsically parallel 
and and the data structure aspects can be explained more succinctly via an explicit 
algorithm. It should, however, be mentioned that the choice of the solution algorithm 
can greatly influence the design and implementation of the necessary mapping and data 
structure. 

Consider the explicit integration of the equations of motion for the structure (2a) 
as recalled here: 

'Mq T f ml = f T f cont (3) 

where f m( and f con < are the internal and applied control forces, respectively, given by 

f int = Dq 4- (K, 4 K a )q 

feont = S® 

The use of the central difference algorithm to integrate (3) leads to the following dif- 
ference equations in time 


q" + * = q"-* + hM -1 (f" 4 t n eont - f? nt ) 
q n+ 1 = q" 4 hq n+ » 


(4) 
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(Processor Configuration) 


Fig. 3 Physical Domain and Its Mapping Onto 
Hypercube Processors 
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4. Implementation and Illustrative Example 


The mapping, partitioning and data structures above discussed have been imple- 
mented based on a shared-memory concurrent machine (Alliant FX/8) by modifying 
the software framework developed for finite element computations [26] and the control- 
structure interaction simulation and design software developed in [27, 28]. At present 
the following specialized systems of equations are implemented: 


where 


and 


Structure: a) Mq + Dq + Kq = f + Bu + Gw 

q(0) = q<» q(o) = qo 

Sensor Output: b) z = Hx +■ m 

Estimator: c) x = Ax + Ef + Bu + L(z - Hx) 

x(0) = 0 

l Control Force: d ) u = — Fx 

H = [H„ H e ] , L =[lj]’ r =l F i F *] 


( 6 ) 


It is noted that in the above implemented equations, we have merged the actua- 
tor and the control law equations into one by neglecting the actuator and control law 
dynamics. Instead, we have introduced an estimator equation as we do not have all 
the measurements needed for complete feedback. In the above equations, B and ,3 
represent the input influence matrix for actuator locations whereas G and G represent 
the disturbance locations. The vector q is the generalized displacement, w is a distur- 
bance vector and the vector m is measurement noise. In Eq. (6b), z is the measured 
sensor output. The matrix H^ is the matrix of displacement sensor locations and H„ 
is the matrix of velocity sensor locations. The state estimator in Eq. (6c) may or may 
not be model based. The superscript ' and ' denote the estimated states and time 
differentiation respectively. The input command, u, is a function of the state estimator 
variables, q and q, and Fi and Fj are control gains. The observer is governed by A, 
the state matrix representing the plant dynamics, and L, the filter gain matrix. 

The software thus implemented was used to test its applicability to solve the 
control-structure interaction design of a model Earth Pointing Satellite (EPS), shown 
in Fig. 6, which is a derivative of a geostationary platform proposed for the study of 
Earth Observation Sciences. Two flexible antennas are attached to a truss bus. Typ- 
ical missions involve pointing one antenna to earth, while tracking or scanning with 



Fig. G 


Earth Pointing Satellite Structure 
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Table 1. EPS Vibration Frequencies (Hz.) 
Mode No. Frequency 


(1-6) 

. 0.000 

( 7 ) 

0.242 

(8) 

0.406 

( 9 ) 

0.565 

(10) 

0.656 

(11,12) 

0.888 

(13) 

1.438 

(14) 

1.536 

(15,16) 

1.776 

(17,18) 

3.026 

(19) 

3.513 

(20) 

3.531 


A small disturbance force was applied to the nominal EPS system in the form of a 
reboost maneuver. The force acted at the center of gravity in the Y-axis direction for 
0 seconds at a 10 N force level and from 0.1 to 0.2 seconds the force level was -10 N. 
The disturbance was removed after 0.2 seconds. Figure 7 shows the open-loop angular 
response about the X-axis of the 15 m antenna. A small amount of passive damping was 
assumed (D = 0.0002 K). The vibrational response produced more than 4.5 n radians 
of RMS pointing error due to this small reboost disturbance. Although many modes 
participate in the flexible body response, this particular reboost maneuver strongly 
excites modes near 4 Hz. The following paragraphs present an integrated control and 
structure design which seeks to lower the vibrational response of the EPS subject to 
some additional constraints. Figure 8 shows the closed-loop angular response about the 
X-axis of the 15 m antenna after design optimization. The pointing error is significantly 
reduced from that of the open-loop system shown . 


5. Future Work and Discussions 

The example problem analyzed in the previous section used a set of lumped actuators 
and localized sensors instead of distributed adaptive actuators and spatially integrated 
sensors. While such a model at best capture the adpative elements used by Anderson 
et al. [29], Matsunaga [30], and Takahara [31], it can not simulate on a large scale 
the distributed usage of piezoelectric actuators and sensors proposed by de Luis [32], 
Rogers et al. [33], and Burk and Hubbard [34], Our immediate future work will 
concentrate on the implementation of distributed adpative elements and assess their 
practical applicability beyond the currently reported beam-like structural components. 
In this regard, we are exploring an adaptation of neural-network concepts [35] in the 
modeling and parallel computations of controlled structures with adaptive elements. 
Specifically, the limits of the applicability of distributed parameter modeling and control 
theory and discrete structures with discrete actuators and sensors, and their cross-over 



Fig. 7 Open-Loop Response of EPS Structure 


x I0' s EPS Antenna Pointing: Closed-Loop Response 



Fig. 8 Closed-Loop Response of Structure EPS 
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performance must be investigated. Design, modeling, simulation and testing criteria 
from such studies will provide greater insight into the eventual adoptions of adaptive 
structures as viable choice for future space systems design alternatives. 

The real-time simulation procedures presented herein may be applicable to the vi- 
bration control of lifeline equipment, and secondarily in minimizing the damage of build- 
ings during earthquakes. In this applications, the sensor measurements used herein can 
be directly applicable to the vibration and earthquake-causing forces on the structures. 
An idea that may prove to be crucial in this case is the use of earthquake-generated 
natural force as vibration minimization actuatori forces. In other words, instead of 
trying to mitigate the earthquake-generating forces, exploit the natural forces instantly 
to activate certain vibration minimizing devices! Research along this line may in the 
end lead to the design of actuators attachable to the columns and floors, if properly 
triggered during earthquakes, can minimize damages based on the natural forces. 
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